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
56IJK_Interfaces::IJK_Interfaces()
65 int size_fixed_arr[3] = {0,0,0};
66 for (
int c=0; c<3; c++)
67 size_fixed_arr[c] = fixed_arr[c].size_array();
68 assert(size_fixed_arr[0] == size_fixed_arr[1]);
69 assert(size_fixed_arr[0] == size_fixed_arr[2]);
71 tab.
resize(size_fixed_arr[0],3);
72 for (
int c=0; c<3; c++)
73 for (
int l=0; l<size_fixed_arr[0]; l++)
74 tab(l,c) = fixed_arr[c](l);
95 Nom fdsom = prefix +
Nom(
"SOMMETS");
96 Nom fdelem = prefix +
Nom(
"ELEMENTS");
97 Nom fdpelocal = prefix +
Nom(
"FACETTE_PE_LOCAL");
98 Nom fdpeowner = prefix +
Nom(
"FACETTE_PE_OWNER");
99 Nom fd_sommet_reel = prefix +
Nom(
"INDEX_SOMMET_REEL");
104 int nsomtot = nb_som;
105 int nelemtot = nb_elem;
112 FloatTab tmp(nb_som, 3);
113 const DoubleTab& sommets = mesh.
sommets();
115 for (
int i = 0; i < nb_som; i++)
116 for (
int j = 0; j < 3; j++)
117 tmp(i, j) = (float)sommets(i, j);
131 IntTab tmp_facettes_renum(nb_elem, 3);
132 const IntTab& facettes = mesh.
facettes();
133 for (
int i = 0; i < nb_elem; i++)
134 for (
int j = 0; j < 3; j++)
135 tmp_facettes_renum(i, j) = facettes(i, j) + offset_sommet;
143 ArrOfInt tmp2(nb_elem);
145 file.
put(tmp2.
addr(), nb_elem, 1);
152 file.
put(tmp2.
addr(), nb_elem, 1);
159 ArrOfInt indice_global(nb_som);
162 ArrOfInt offset_received(np);
163 ArrOfInt offset_to_send(np);
166 offset_to_send = offset_sommet;
168 envoyer_all_to_all(offset_to_send, offset_received);
175 for (
int i = 0; i < nb_som; i++)
177 const int indice_local_sur_proc_proprietaire = sommet_num_owner[i];
178 const int proc_proprietaire = sommet_PE_owner[i];
179 const int offset_proc_proprietaire = offset_received[proc_proprietaire];
180 indice_global[i] = indice_local_sur_proc_proprietaire + offset_proc_proprietaire;
182 file.
ouvrir(fd_sommet_reel);
188 const Nom format =
"INT32";
192 master_file.
ouvrir(filename, ios::app);
193 master_file <<
"Geom " << meshname <<
" type_elem=TRIANGLE_3D" << finl;
195 master_file <<
"Champ SOMMETS " << basename(fdsom) <<
" geometrie=" << meshname;
196 master_file <<
" size=" << nsomtot <<
" composantes=3" << finl;
198 master_file <<
"Champ ELEMENTS " << basename(fdelem) <<
" geometrie=" << meshname;
199 master_file <<
" size=" << nelemtot <<
" composantes=3 format=" << format << finl;
201 master_file <<
"Champ FACETTE_PE_LOCAL " << basename(fdpelocal) <<
" geometrie=" << meshname;
202 master_file <<
" size=" << nelemtot <<
" composantes=1 format=" << format <<
" localisation=ELEM" << finl;
204 master_file <<
"Champ FACETTE_PE_OWNER " << basename(fdpeowner) <<
" geometrie=" << meshname;
205 master_file <<
" size=" << nelemtot <<
" composantes=1 format=" << format <<
" localisation=ELEM" << finl;
207 master_file <<
"Champ INDEX_SOMMET_REEL " << basename(fd_sommet_reel) <<
" geometrie=" << meshname;
208 master_file <<
" size=" << nsomtot <<
" composantes=1 format=" << format <<
" localisation=SOM" << finl;
213void runge_kutta3_update(
const DoubleTab& dvi, DoubleTab& G, DoubleTab& l,
const int step,
const double dt_tot,
216 const double coeff_a[3] = {0., -5. / 9., -153. / 128.};
218 const double coeff_Gk[3] = {1., 4. / 9., 15. / 32.};
220 const double facteurG = coeff_a[step];
221 const double intermediate_dt = compute_fractionnal_timestep_rk3(dt_tot, step);
222 const double delta_t_divided_by_Gk = intermediate_dt / coeff_Gk[step];
234 for (
int isom = 0; isom < nbsom; isom++)
236 for (
int dir = 0; dir < 3; dir++)
240 G(isom, dir) = 111. / intermediate_dt;
241 l(isom, dir) = 111. / intermediate_dt;
243 double x = dvi(isom, dir);
245 l(isom, dir) += x * delta_t_divided_by_Gk;
251 for (
int isom = 0; isom < nbsom; isom++)
253 for (
int dir = 0; dir < 3; dir++)
257 G(isom, dir) = 111. / intermediate_dt;
258 l(isom, dir) = 111. / intermediate_dt;
260 double x = G(isom, dir) * facteurG + dvi(isom, dir);
262 l(isom, dir) += x * delta_t_divided_by_Gk;
268 for (
int isom = 0; isom < nbsom; isom++)
270 for (
int dir = 0; dir < 3; dir++)
274 G(isom, dir) = 111. / intermediate_dt;
275 l(isom, dir) = 111. / intermediate_dt;
277 double x = G(isom, dir) * facteurG + dvi(isom, dir);
278 l(isom, dir) += x * delta_t_divided_by_Gk;
283 Cerr <<
"Error in runge_kutta_update: wrong step" << finl;
297 os <<
" terme_gravite rho_g \n";
299 os <<
" terme_gravite grad_i \n";
302 os <<
" correction_gradient_potentiel \n";
314 os <<
" active_repulsion_paroi \n";
320 Cerr <<
"IJK_Interfaces::printOn -- couleurs des bulles stockees " << finl;
321 os <<
" follow_colors " <<
"\n"
326 Cerr <<
"IJK_Interfaces::printOn -- Groupes d'appartenance stockes " << finl;
332 if (
parcours_.get_correction_parcours_thomas())
333 os <<
" parcours_interface { correction_parcours_thomas } " <<
"\n";
334 if (
parcours_.get_parcours_sans_tolerance())
335 os <<
" parcours_interface { parcours_sans_tolerance } " <<
"\n";
338 os <<
"use_barycentres_velocity" <<
"\n";
340 os <<
"use_barycentres_velocity" <<
"\n";
342 double max_force_compo = 0.;
345 if (max_force_compo > 1.e-16)
348 Cerr <<
"IJK_Interfaces::printOn -- Storing mean into sauv for future "
349 "restart. (max_force_compo= "
350 << max_force_compo <<
" )." << finl;
394 param.ajouter_flag(
"frozen", &
frozen_);
397 param.ajouter(
"parcours_interface",&
parcours_);
414 param.lire_avec_accolades(is);
416 Cout <<
"IJK_Interfaces::readOn : Option gravite : " <<
terme_gravite_
428 Cout <<
"IJK_Interfaces::readOn : il y a : " <<
nb_groups_
429 <<
" classes/groupes de bulles suivis " << finl;
435 <<
" bulles reelles reprises avec une position de reference. " << finl;
439 Cout <<
"IJK_Interfaces::readOn : " <<
mean_force_.dimension(0)
440 <<
" bulles reelles reprises avec une force_moyenne imposee. " << finl;
447 if (un_mot==
"parcours_interface")
454 Cerr <<
"Unknown Keyword " << un_mot <<
" in IJK_Interfaces::lire_motcle_non_standard" << finl;
473 if (!
parcours_.get_correction_parcours_thomas())
475 Cerr <<
"Warning: The cut-cell method overrides default or user parameters to force the use of correction_parcours_thomas in parcours_interface." << finl;
476 parcours_.set_correction_parcours_thomas();
479 if (!
parcours_.get_parcours_sans_tolerance())
481 Cerr <<
"Warning: The cut-cell method overrides default or user parameters to force the use of parcours_sans_tolerance in parcours_interface." << finl;
498 double min_indicatrice = 1.;
501 for (
int n = 0; n < cut_cell_disc.
get_n_loc(); n++)
503 Int3 ijk = cut_cell_disc.
get_ijk(n);
511 if (indicatrice > 0.)
512 min_indicatrice = std::min(min_indicatrice, indicatrice);
514 if (indicatrice <= .5 && indicatrice > .4)
516 else if (indicatrice <= .4 && indicatrice > .3)
518 else if (indicatrice <= .3 && indicatrice > .2)
520 else if (indicatrice <= .2 && indicatrice > .1)
522 else if (indicatrice <= .1 && indicatrice > .05)
524 else if (indicatrice <= .05 && indicatrice > .01)
526 else if (indicatrice <= .01 && indicatrice > .0)
530 assert(indicatrice == 0.);
535 mp_sum_for_each(count_total, count_40pct, count_30pct, count_20pct, count_10pct);
549 Cerr <<
"Bilan de l'indicatrice:" << finl;
550 Cerr <<
" Valeur minimale: " << min_indicatrice << finl;
551 Cerr <<
" Compte sur les cellules de la structure diphasique:" << finl;
552 Cerr <<
" 40-50% " << count_40pct << finl;
553 Cerr <<
" 30-40% " << count_30pct << finl;
554 Cerr <<
" 20-30% " << count_20pct << finl;
555 Cerr <<
" 10-20% " << count_10pct << finl;
556 Cerr <<
" 5-10% " << count_5pct << finl;
557 Cerr <<
" 1-5% " << count_1pct << finl;
558 Cerr <<
" 0-1% " << count_0pct << finl;
559 Cerr <<
" pure " << count_pure << finl;
560 Cerr <<
" ------ " << finl;
561 Cerr <<
" total " << count_total << finl;
569 int iteration_solver_surface_efficace_face = 0;
571 if (type_surface_efficace_face == TYPE_SURFACE_EFFICACE_FACE::CONSERVATION_VOLUME_ITERATIF)
578 iteration_solver_surface_efficace_face,
590 iteration_solver_surface_efficace_face,
605 (type_surface_efficace_interface == TYPE_SURFACE_EFFICACE_INTERFACE::EXPLICITE),
625 if (type_surface_efficace_interface == TYPE_SURFACE_EFFICACE_INTERFACE::CONSERVATION_VOLUME)
652 (type_surface_efficace_face == TYPE_SURFACE_EFFICACE_FACE::EXPLICITE),
662 (type_surface_efficace_interface == TYPE_SURFACE_EFFICACE_INTERFACE::EXPLICITE),
676 const IJK_Field_vector3_double& velocity_ft = ref_ijk_ft_->eq_ns().get_velocity_ft();
678 const DoubleTab& sommets = mesh.
sommets();
680 ArrOfDouble vinterp_component(nbsom);
682 for (
int direction = 0; direction < 3; direction++)
688 ijk_interpolate_skip_unknown_points(velocity_ft[direction], sommets, vinterp_component,
690 for (
int i = 0; i < nbsom; i++)
691 vinterp_(i, direction) = vinterp_component[i];
713 if (nom==
"CONCENTRATION_INTERFACE")
718 if (nom==
"GRADX_CONCENTRATION_INTERFACE")
720 IJK_Field_double& dconc_interf_x =
scalar_post_fields_.at(
"GRADX_CONCENTRATION_INTERFACE");
723 if (nom==
"GRADY_CONCENTRATION_INTERFACE")
725 IJK_Field_double& dconc_interf_y =
scalar_post_fields_.at(
"GRADY_CONCENTRATION_INTERFACE");
728 if (nom==
"GRADZ_CONCENTRATION_INTERFACE")
730 IJK_Field_double& dconc_interf_z =
scalar_post_fields_.at(
"GRADZ_CONCENTRATION_INTERFACE");
733 if (nom==
"LAPLACIAN_CONCENTRATION_INTERFACE")
735 IJK_Field_double& lapla_interf =
scalar_post_fields_.at(
"LAPLACIAN_CONCENTRATION_INTERFACE");
743 if (nom==
"GRADX_SIGMA")
746 dsigma_x.
data() =
maillage_ft_ijk_.Surfactant_facettes().get_interfacial_source_term_sommets(0);
748 if (nom==
"GRADY_SIGMA")
751 dsigma_y.
data() =
maillage_ft_ijk_.Surfactant_facettes().get_interfacial_source_term_sommets(1);
753 if (nom==
"GRADZ_SIGMA")
756 dsigma_z.
data() =
maillage_ft_ijk_.Surfactant_facettes().get_interfacial_source_term_sommets(2);
758 if (nom==
"DISTANCE_AUTRES_INTERFACES")
760 DoubleTab vr_to_closer;
766 return champs_compris_.get_champ(nom);
768 Cerr <<
"ERROR in IJK_Interfaces::get_IJK_field_vector : " << finl;
769 Cerr <<
"Requested field '" << nom <<
"' is not recognized by IJK_Interfaces::get_IJK_field_vector()." << finl;
776 return champs_compris_.get_champ_vectoriel(nom);
778 Cerr <<
"ERROR in IJK_Interfaces::get_IJK_field_vector : " << finl;
779 Cerr <<
"Requested field '" << nom <<
"' is not recognized by IJK_Interfaces::get_IJK_field_vector()." << finl;
785 std::vector<FieldInfo_t> c =
788 {
"INDICATRICE", Entity::ELEMENT, Nature_du_champ::scalaire,
false },
789 {
"INDICATRICE_FT", Entity::ELEMENT, Nature_du_champ::scalaire,
false },
790 {
"OLD_INDICATRICE", Entity::ELEMENT, Nature_du_champ::scalaire,
false },
791 {
"OLD_INDICATRICE_FT", Entity::ELEMENT, Nature_du_champ::scalaire,
false },
792 {
"REPULSION_FT", Entity::ELEMENT, Nature_du_champ::scalaire,
false },
793 {
"CUT_FIELDS_BARY_L", Entity::ELEMENT, Nature_du_champ::vectoriel,
false },
795 {
"COURBURE", Entity::NODE, Nature_du_champ::scalaire,
true },
796 {
"CONCENTRATION_INTERFACE", Entity::ELEMENT, Nature_du_champ::scalaire,
true },
797 {
"GRADX_CONCENTRATION_INTERFACE", Entity::NODE, Nature_du_champ::scalaire,
true },
798 {
"GRADY_CONCENTRATION_INTERFACE", Entity::NODE, Nature_du_champ::scalaire,
true },
799 {
"GRADZ_CONCENTRATION_INTERFACE", Entity::NODE, Nature_du_champ::scalaire,
true },
800 {
"SIGMA", Entity::NODE, Nature_du_champ::scalaire,
true },
801 {
"GRADX_SIGMA", Entity::NODE, Nature_du_champ::scalaire,
true },
802 {
"GRADY_SIGMA", Entity::NODE, Nature_du_champ::scalaire,
true },
803 {
"GRADZ_SIGMA", Entity::NODE, Nature_du_champ::scalaire,
true },
804 {
"LAPLACIAN_CONCENTRATION_INTERFACE", Entity::ELEMENT, Nature_du_champ::scalaire,
true },
805 {
"DISTANCE_AUTRES_INTERFACES", Entity::NODE, Nature_du_champ::scalaire,
true }
808 chps.insert(chps.end(), c.begin(), c.end());
813 for (
const auto& n : champs_compris_.liste_noms_compris())
815 for (
const auto& n : champs_compris_.liste_noms_compris_vectoriel())
821 Cerr <<
"Reset bubble rising velocity calculations" << finl;
823 Cerr <<
"Compute compo_connex from bounding box" << finl;
825 Cerr <<
"Compute compo_connex from interface compo in mixed cells" << finl;
827 Cerr <<
"Compute rising velocity from compo connex (barycentre calc)" << finl;
834 const int thermal_probes_ghost_cells,
835 const bool compute_vint,
836 const bool is_switch)
838 Cerr <<
"Entree dans IJK_Interfaces::initialize" << finl;
841 ref_domaine_ = domaine_FT;
851 const int nb_ghost_cells = std::max(thermal_probes_ghost_cells, (
int) 4);
901 allocate_cell_vector(groups_indicatrice_ft_test_, domaine_FT, 1);
902 allocate_cell_vector(groups_indicatrice_ft_test_, domaine_FT, 1);
906 for (
int i = 0; i < max_authorized_nb_of_components_; i++)
914 for (
int dir = 0; dir < 3; dir++)
922 for (
int dir = 0; dir < 3; dir++)
927 for (
int dir = 0; dir < 3; dir++)
929 const int idx = i * 3 + dir;
954 for (
int bary_compo = 0; bary_compo < 3; bary_compo++)
967 for (
int face_dir = 0; face_dir < 3; face_dir++)
980 for (
int face_dir = 0; face_dir < 3; face_dir++)
988 for (
int face_dir = 0; face_dir < 3; face_dir++)
990 for (
int bary_compo = 0; bary_compo < 2; bary_compo++)
999 for (
int face_dir = 0; face_dir < 3; face_dir++)
1001 for (
int bary_compo = 0; bary_compo < 2; bary_compo++)
1025 for (
int d = 0; d < 3; d++)
1039 Cut_field_vector3_double& cut_field_deformation_velocity =
static_cast<Cut_field_vector3_double&
>(
deformation_velocity_);
1040 allocate_velocity_ephemere(*ref_ijk_ft_->get_cut_cell_disc(), cut_field_deformation_velocity, domaine_NS, 2);
1063 for (
int direction = 0; direction < 3; direction++)
1065 const double ori = domaine_NS.
get_origin(direction);
1075 const double oriFT = domaine_FT.
get_origin(direction);
1083 double duplicate_stencil_width ;
1086 duplicate_stencil_width =
1090 duplicate_stencil_width =
1170 DoubleTab centre_gravite;
1171 ArrOfDouble volume_bulles;
1189 DoubleTab vr_to_closer;
1209 Cerr <<
"Error in option combination: invalid choice of optimized methods "
1210 <<
"for both the color_function in the forces computation "
1211 <<
"(parser_=0) and the indicator function calculation "
1212 <<
"(recompute_indicator_=0)"
1214 Cerr <<
"Maybe you are not using the forces computation?" << finl;
1223 const Domaine& domaine = domaine_vf.
domaine();
1238 courb.nommer(
"COURBURE");
1239 champs_compris_.ajoute_champ(courb);
1243 conc_interf.nommer(
"CONCENTRATION_INTERFACE");
1244 champs_compris_.ajoute_champ(conc_interf);
1248 dconc_interf_x.nommer(
"GRADX_CONCENTRATION_INTERFACE");
1249 champs_compris_.ajoute_champ(dconc_interf_x);
1253 dconc_interf_y.nommer(
"GRADY_CONCENTRATION_INTERFACE");
1254 champs_compris_.ajoute_champ(dconc_interf_y);
1258 dconc_interf_z.nommer(
"GRADZ_CONCENTRATION_INTERFACE");
1259 champs_compris_.ajoute_champ(dconc_interf_z);
1263 sigma.nommer(
"SIGMA");
1264 champs_compris_.ajoute_champ(sigma);
1268 dsigma_x.nommer(
"GRADX_SIGMA");
1269 champs_compris_.ajoute_champ(dsigma_x);
1273 dsigma_y.nommer(
"GRADY_SIGMA");
1274 champs_compris_.ajoute_champ(dsigma_y);
1278 dsigma_z.nommer(
"GRADZ_SIGMA");
1279 champs_compris_.ajoute_champ(dsigma_z);
1283 lapla_interf.nommer(
"LAPLACIAN_CONCENTRATION_INTERFACE");
1284 champs_compris_.ajoute_champ(lapla_interf);
1288 d.nommer(
"DISTANCE_AUTRES_INTERFACES");
1289 champs_compris_.ajoute_champ(d);
1294 Cerr <<
"IJK_Interfaces::milieu not coded ! access it from NS. " << finl;
1300 Cerr <<
"IJK_Interfaces::milieu not coded ! access it from NS. " << finl;
1308 Cerr <<
"Error for the method IJK_Interfaces::associer_pb_base\n";
1309 Cerr <<
" IJK_Interfaces equation must be associated to\n";
1310 Cerr <<
" a Probleme_FTD_IJK_base problem type\n";
1327 ref_ijk_ft_switch_ = ijk_ft_switch;
1336 liste.add(
"INTERFACES");
1337 liste.add(
"COMPO_CONNEXE");
1338 liste.add(
"DISTANCE_AUTRES_INTERFACES");
1340 liste.add(
"COLOR_Y");
1345 liste.add(
"RK3_STORE_VI");
1352 const char *lata_name,
1353 const int lata_step)
const
1358 if (liste_post_instantanes.
contient_(
"INTERFACES"))
1362 if (liste_post_instantanes.
contient_(
"COLOR_Y"))
1366 Cerr <<
"On demande le post-traitement du champ COLOR_Y alors que le flag follow_colors "
1367 <<
"n'a pas ete active... Ce champ est donc inconnu. Post-traitement refuse. " << finl;
1372 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"COLOR_Y",
"ELEM", color, lata_step);
1374 if (liste_post_instantanes.
contient_(
"VINTERP"))
1376 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"VINTERP",
"SOM",
RK3_G_store_vi_, lata_step);
1378 if (liste_post_instantanes.
contient_(
"VI"))
1387 for (
int dir= 0; dir < 3; dir++)
1391 Cerr <<
"Posttraitement du champ VI vide : on stocke 0." << finl;
1395 for (
int i = 0; i < nbsom; i++)
1399 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"VI_X",
"SOM", vi, lata_step);
1401 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"VI_Y",
"SOM", vi, lata_step);
1403 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"VI_Z",
"SOM", vi, lata_step);
1408 if (liste_post_instantanes.
contient_(
"RK3_STORE_VI"))
1410 ArrOfDoubleFT store_vi_compo;
1417 for (
int dir= 0; dir < 3; dir++)
1422 Cerr <<
"Posttraitement du champ RK3_STORE_VI vide : on stocke 0." << finl;
1423 store_vi_compo = 0.;
1427 for (
int i = 0; i < nbsom; i++)
1430 store_vi_compo[i] = val;
1434 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"RK3_STORE_VI_X",
"SOM", store_vi_compo, lata_step);
1436 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"RK3_STORE_VI_Y",
"SOM", store_vi_compo, lata_step);
1438 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"RK3_STORE_VI_Z",
"SOM", store_vi_compo, lata_step);
1459 const Domaine_IJK& geom_FT = ref_domaine_.valeur();
1460 for (
int direction = 0; direction < 3; direction++)
1466 const double oriFT = geom_FT.
get_origin(direction);
1484 DoubleTab bounding_box;
1490 ArrOfInt masque_delete_pour_compo;
1500 int nb_bulles_reelles_deleted = 0;
1501 ArrOfInt old_to_new_compo;
1506 for (
int i = 0; i < masque_delete_pour_compo.
size_array(); i++)
1508 if (masque_delete_pour_compo[i] != 0)
1511 nb_bulles_reelles_futur -= 1;
1512 nb_bulles_reelles_deleted += 1;
1513 old_to_new_compo[i] = -1;
1517 old_to_new_compo[i] = count;
1521 nb_bulles_reelles_futur = count;
1524 envoyer_broadcast(old_to_new_compo, 0);
1525 envoyer_broadcast(nb_bulles_reelles_futur, 0);
1526 envoyer_broadcast(flag, 0);
1527 envoyer_broadcast(nb_bulles_reelles_deleted, 0);
1531 Cerr << nb_bulles_reelles_deleted <<
" marked for deletion out of " <<
nb_bulles_reelles_
1532 <<
" so that only " << nb_bulles_reelles_futur <<
" will remain." << finl;
1536 for (
int fa7 = 0; fa7 < mesh.
nb_facettes(); fa7++)
1538 icompo = compo_connexe_facettes[fa7];
1539 assert(icompo >= 0);
1550 <<
" to " << nb_bulles_reelles_futur << finl;
1558 int inew = old_to_new_compo[icompo];
1559 assert(inew <= icompo);
1576 Cerr <<
"The table of bubbles groups (compo_to_group_), "
1577 <<
" (as well as possibly through_yminus_ and positions_reference_ "
1578 <<
"if needed)" <<
" have been updated in accordance." << finl;
1611 os <<
"# Impression des couleurs de bulles" << finl;
1612 os <<
"# colonne 1 : temps" << finl;
1613 os <<
"# colonne K : Couleur de la bulle K-2" << finl;
1615 os << current_time <<
" ";
1616 for (
int j = 0; j < n; j++)
1627 for (
int i = 0; i < n; i++)
1633 const int compo = compo_facettes[i];
1656 for (
int i = 0; i < n; i++)
1664 Cerr <<
"Exception dans IJK_Interfaces::get_ghost_number_from_compo(). "
1665 <<
" Compo demandee " << compo <<
" introuvable...";
1676 surfaces.
resize_array(nbulles_reelles + nbulles_ghost, RESIZE_OPTIONS::NOCOPY_NOINIT);
1680 for (
int i = 0; i < n; i++)
1684 int compo = compo_facettes[i];
1691 compo = nbulles_reelles - 1 - idx_ghost;
1694 const double s = surfaces_facettes[i];
1695 surfaces[compo] += s;
1706 ArrOfDouble& out)
const
1714 for (
int i = 0; i < n; i++)
1718 int compo = compo_facettes[i];
1723 const double s = in[i];
1728 for (
int c = 0; c < nbulles; c++)
1729 out[c] /= surfaces[c];
1739 Cerr <<
"interfacial surface and normale are sure to be up-to-date" << finl;
1757 bool tmp_readen_flag =
false;
1758 Nom suffix_tmp =
".old";
1761 if (tmp_readen_flag)
1765 suffix_tmp =
".new";
1768 if (tmp_readen_flag)
1775 if (tmp_readen_flag)
1787 ArrOfDouble& bubbles_rising_vel_mag)
1789 bool is_readen =
false;
1792 bubbles_rising_vel_mag.
reset();
1793 for (
int c=0; c<3; c++)
1795 bubbles_rising_dir[c].reset();
1796 bubbles_rising_vel[c].reset();
1798 if (bubbles_bary_computed)
1800 int line_counter = 0;
1804 const Nom interf_dir = dirname(interf_name);
1806 const Nom suffix_bary =
".sauv.barycentres";
1807 const Nom file_folder = interf_dir + case_name + suffix_bary +
".vel";
1810 read_tmp.open(file_folder);
1811 const bool read_file = read_tmp ? true :
false;
1815 EFichier fic_bary_vel(file_folder);
1816 ifstream& ifstream_bary_vel = fic_bary_vel.
get_ifstream();
1817 const char delimiter =
' ';
1818 Cerr <<
"Read coordinates of bubbles barycentres from: " << file_folder << finl;
1819 while (std::getline(ifstream_bary_vel, line))
1821 std::stringstream ssline(line);
1822 Cerr <<
"Line number: " << line_counter << finl;
1823 while (std::getline(ssline, line, delimiter))
1827 Cerr <<
"Param: " << line << finl;
1835 dir = var_index - 1;
1836 bubbles_rising_vel[dir].append_array(std::stod(line));
1841 dir = var_index - 4;
1842 bubbles_rising_dir[dir].append_array(std::stod(line));
1856 fic_bary_vel.
close();
1865 bool is_readen =
false;
1868 for (
int c=0; c<3; c++)
1870 bubbles_bary[c].reset();
1872 if (bubbles_bary_computed)
1874 int line_counter = 0;
1878 const Nom interf_dir = dirname(interf_name);
1880 const Nom suffix_bary =
".sauv.barycentres";
1881 const Nom file_folder = interf_dir + case_name + suffix_bary + suffix;
1883 read_tmp.open(file_folder);
1884 const bool read_file = read_tmp ? true :
false;
1890 const char delimiter =
' ';
1891 Cerr <<
"Read coordinates of bubbles barycentres from: " << file_folder << finl;
1892 while (std::getline(ifstream_bary_old, line))
1894 std::stringstream ssline(line);
1895 Cerr <<
"Line number: " << line_counter << finl;
1896 while (std::getline(ssline, line, delimiter))
1900 Cerr <<
"Param: " << line << finl;
1908 dir = var_index - 1;
1909 bubbles_bary[dir].append_array(std::stod(line));
1929 const Nom end_space =
" ";
1930 const Nom escape =
"\n";
1935 const Nom interf_dir = dirname(interf_name);
1936 const Nom suffix =
".sauv.barycentres";
1937 const int reset = 1;
1938 const Nom bary_header_old =
"ibubble bary_old_x bary_old_y bary_old_z";
1939 const Nom bary_header_new =
"ibubble bary_new_x bary_new_y bary_new_z";
1940 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";
1942 SFichier fic_bary_old = Open_file_folder(interf_dir, suffix +
".old", bary_header_old, reset, 0);
1943 for (
int ibubble=0; ibubble<nbulles_reelles; ibubble++)
1945 fic_bary_old << ibubble << end_space;
1950 fic_bary_old.
close();
1951 SFichier fic_bary_new = Open_file_folder(interf_dir, suffix +
".new", bary_header_new, reset, 0);
1952 for (
int ibubble=0; ibubble<nbulles_reelles; ibubble++)
1954 fic_bary_new << ibubble << end_space;
1959 fic_bary_new.
close();
1960 SFichier fic_bary_vel = Open_file_folder(interf_dir, suffix +
".vel", bary_vel_header, reset, 0);
1961 for (
int ibubble=0; ibubble<nbulles_reelles; ibubble++)
1963 fic_bary_vel << ibubble << end_space;
1972 fic_bary_vel.
close();
1977 DoubleTab& barycentres,
1978 const int& store_values)
1984 if (ref_ijk_ft_->schema_temps_ijk().get_tstep() == 0)
1997 for (
int i = 0; i < nbulles_reelles; i++)
2000 for (
int dir=0; dir<3; dir++)
2005 const int old_new_up_down = (pos_new - pos_old) < (-ldir / 2.);
2006 const int old_new_down_up = (pos_new - pos_old) > (ldir / 2.);
2007 if (old_new_up_down || old_new_down_up)
2009 if (old_new_up_down)
2022 for (
int dir=0; dir<3; dir++)
2035 DoubleTab& centre_gravite)
const
2041 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2042 volumes.
resize_array(nbulles_tot, RESIZE_OPTIONS::NOCOPY_NOINIT);
2044 centre_gravite.
resize(nbulles_tot, 3, RESIZE_OPTIONS::NOCOPY_NOINIT);
2045 centre_gravite = 0.;
2048 const IntTab& facettes = mesh.
facettes();
2049 const DoubleTab& sommets = mesh.
sommets();
2051 for (
int i = 0; i < n; i++)
2055 int compo = compo_facettes[i];
2062 compo = nbulles_reelles - 1 - idx_ghost;
2064 const double s = surfaces_facettes[i];
2065 const double normale_scalaire_direction = normales_facettes(i, 0);
2067 const int i0 = facettes(i, 0);
2068 const int i1 = facettes(i, 1);
2069 const int i2 = facettes(i, 2);
2070 const double coord_centre_gravite_i = (sommets(i0, 0) + sommets(i1, 0) + sommets(i2, 0)) / 3.;
2071 const double coord_centre_gravite_j = (sommets(i0, 1) + sommets(i1, 1) + sommets(i2, 1)) / 3.;
2072 const double coord_centre_gravite_k = (sommets(i0, 2) + sommets(i1, 2) + sommets(i2, 2)) / 3.;
2073 const double volume_prisme = coord_centre_gravite_i * s * normale_scalaire_direction;
2075 centre_gravite(compo, 0) += volume_prisme * (coord_centre_gravite_i * 0.5);
2076 centre_gravite(compo, 1) += volume_prisme * coord_centre_gravite_j;
2077 centre_gravite(compo, 2) += volume_prisme * coord_centre_gravite_k;
2078 volumes[compo] += volume_prisme;
2083 for (
int i = 0; i < nbulles_tot; i++)
2086 const double x = (volumes[i] == 0.) ? 0. : 1. / volumes[i];
2087 centre_gravite(i, 0) *= x;
2088 centre_gravite(i, 1) *= x;
2089 centre_gravite(i, 2) *= x;
2100 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2101 const IntTab& facettes = mesh.
facettes();
2102 const DoubleTab& sommets = mesh.
sommets();
2106 ArrOfDouble volumes;
2107 DoubleTab centre_gravite;
2110 DoubleTab d_max(nbulles_tot);
2112 DoubleTab d_min(nbulles_tot);
2118 for (
int i = 0; i < n; i++)
2122 int compo = compo_facettes[i];
2129 compo = nbulles_reelles - 1 - idx_ghost;
2132 const int i0 = facettes(i, 0);
2133 const int i1 = facettes(i, 1);
2134 const int i2 = facettes(i, 2);
2135 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)) );
2136 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)) );
2137 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)) );
2140 d_imax = std::max(d_i_0,std::max(d_i_1,d_i_2));
2141 d_imin = std::min(d_i_0,std::min(d_i_1,d_i_2));
2144 if (d_imax > d_max[compo])
2145 d_max[compo] = d_imax;
2147 if (d_imin < d_min[compo])
2148 d_min[compo] = d_imin;
2153 for (
int i = 0; i < nbulles_tot; i++)
2155 aspect_ratio[i] = d_max[i]/d_min[i];
2174 for (
int bulle = 0; bulle < nbulles_reelles; bulle++)
2176 surfactant_min(bulle) = 1.e10;
2177 surfactant_max(bulle) = -1.e10;
2179 for (
int i = 0; i < n; i++)
2183 surfactant(compo_facettes(i))+=Surf(i)*Sfa7(i);
2184 if(Surf(i)>surfactant_max(compo_facettes(i)))
2186 surfactant_max(compo_facettes(i)) = Surf(i);
2188 if(Surf(i)<surfactant_min(compo_facettes(i)))
2190 surfactant_min(compo_facettes(i)) = Surf(i);
2202 DoubleTab& poussee)
const
2208 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2209 poussee.
resize(nbulles_tot, 3, RESIZE_OPTIONS::NOCOPY_NOINIT);
2213 const IntTab& facettes = mesh.
facettes();
2214 const DoubleTab& sommets = mesh.
sommets();
2216 for (
int i = 0; i < n; i++)
2220 int compo = compo_facettes[i];
2227 compo = nbulles_reelles - 1 - idx_ghost;
2229 const double s = surfaces_facettes[i];
2231 const int i0 = facettes(i, 0);
2232 const int i1 = facettes(i, 1);
2233 const int i2 = facettes(i, 2);
2234 const double coord_centre_gravite_i = (sommets(i0, 0) + sommets(i1, 0) + sommets(i2, 0)) / 3.;
2235 const double coord_centre_gravite_j = (sommets(i0, 1) + sommets(i1, 1) + sommets(i2, 1)) / 3.;
2236 const double coord_centre_gravite_k = (sommets(i0, 2) + sommets(i1, 2) + sommets(i2, 2)) / 3.;
2237 const double grav_scalaire_position_fois_s =
2238 (grav(0,0) * coord_centre_gravite_i + grav(0,1) * coord_centre_gravite_j + grav(0,2) * coord_centre_gravite_k) * s;
2239 for (
int dir = 0; dir < 3; dir++)
2240 poussee(compo, dir) += grav_scalaire_position_fois_s * normales_facettes(i, dir);
2257 const double vol = dxi * dxj * dxk;
2262 for (
int fa7 = 0; fa7 < n; fa7++)
2269 const double sf = surface_facettes[fa7];
2280 ai(ijk[0], ijk[1], ijk[2]) += surf / vol;
2299 const double vol = dxi * dxj * dxk;
2304 for (
int fa7 = 0; fa7 < n; fa7++)
2310 int compo_fa7 = compo_connex[fa7];
2317 if (compo_fa7 != compo)
2321 const double sf = surface_facettes[fa7];
2332 ai(ijk[0], ijk[1], ijk[2]) += surf / vol;
2351 const double vol = dxi * dxj * dxk;
2353 for (
int fa7 = 0; fa7 < n; fa7++)
2360 const double sf = surface_facettes[fa7];
2361 const int compo_fa7 = compo_connex[fa7];
2363 if (compo_fa7 != compo)
2376 if (ijk[0]==i_ref and ijk[1]==j_ref and ijk[2]==k_ref)
2388static int encoder_compo(
int num_bulle,
int code_deplacement)
2394 return (num_bulle << 6) | (code_deplacement);
2397static int decoder_numero_bulle(
const int code)
2399 const int num_bulle = code >> 6;
2410 const IntTab& facettes = mesh.
facettes();
2411 const ArrOfDouble& courbure =
maillage_ft_ijk_.get_update_courbure_sommets();
2417 IJK_Field_double SI_ft;
2421 kappa_ft.
data() = 0.;
2423 for (
int fa7 = 0; fa7 < n; fa7++)
2429 const double sf=surface_facettes[fa7];
2437 for (
int isom = 0; isom< 3; isom++)
2439 const int num_som = facettes(fa7, isom);
2440 const double kappa = courbure[num_som];
2441 kappa_ft(ijk[0],ijk[1],ijk[2]) += kappa*surf/3.;
2443 SI_ft(ijk[0],ijk[1],ijk[2]) += surf;
2449 const int nx = kappa_ft.
ni();
2450 const int ny = kappa_ft.
nj();
2451 const int nz = kappa_ft.
nk();
2453 for (
int k = 0; k < nz; k++)
2454 for (
int j = 0; j < ny; j++)
2455 for (
int i = 0; i < nx; i++)
2457 if (kappa_ft(i,j,k)!=0.)
2458 kappa_ft(i,j,k)/=SI_ft(i,j,k);
2475 IJK_Field_vector3_double& normale_cell,
2476 const int igroup)
const
2482 const IntTab& facettes = mesh.
facettes();
2483 const ArrOfDouble& courbure =
maillage_ft_ijk_.get_update_courbure_sommets();
2490 const double vol = dxi * dxj * dxk;
2494 kappa_ai.
data() = 0.;
2495 for (
int dir = 0; dir < 3; dir++)
2496 normale_cell[dir].data() = 0.;
2499 for (
int fa7 = 0; fa7 < n; fa7++)
2508 int icompo = compo_facette[fa7];
2512 icompo = decoder_numero_bulle(-icompo);
2517 const double sf=surface_facettes[fa7];
2526 for (
int dir = 0; dir< 3; dir++)
2528 const double nx = normale_facettes(fa7,dir);
2529 normale_cell[dir](ijk[0],ijk[1],ijk[2]) += nx * surf;
2531 const double fac = surf/(3.*vol);
2532 for (
int isom = 0; isom< 3; isom++)
2534 const int num_som = facettes(fa7, isom);
2535 const double kappa = courbure[num_som];
2536 kappa_ai(ijk[0],ijk[1],ijk[2]) += kappa*fac;
2539 ai(ijk[0],ijk[1],ijk[2]) += surf/vol;
2547 const int nx = ai.
ni();
2548 const int ny = ai.
nj();
2549 const int nz = ai.
nk();
2551 for (
int k = 0; k < nz; k++)
2552 for (
int j = 0; j < ny; j++)
2553 for (
int i = 0; i < nx; i++)
2555 double norme_carre = 0.;
2556 for (
int dir = 0; dir < 3; dir++)
2558 const double dnx = normale_cell[dir](i, j, k);
2559 norme_carre += dnx * dnx;
2562 if (norme_carre > 0.)
2563 for (
int dir = 0; dir < 3; dir++)
2564 normale_cell[dir](i, j, k) *= 1. / sqrt(norme_carre);
2571 const ArrOfDouble& surface_facette,
2572 const ArrOfDouble& surface_par_bulle,
2573 const ArrOfInt& compo_connexes_facettes,
2574 const int nbulles_reelles,
2575 const int nbulles_ghost,
2576 const DoubleTab& vitesse_sommets,
2577 DoubleTab& vitesses_translation_bulles)
const
2579 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2580 assert(vitesses_translation_bulles.
dimension(0) == nbulles_tot);
2581 assert(surface_par_bulle.
size_array() == nbulles_tot);
2582 assert(vitesses_translation_bulles.
dimension(1) == 3);
2584 const IntTab& facettes = maillage.
facettes();
2586 vitesses_translation_bulles = 0.;
2594 for (
int fa7 = 0; fa7 < maillage.
nb_facettes(); fa7++)
2601 int compo = compo_connexes_facettes[fa7];
2609 compo = nbulles_reelles - 1 - idx_ghost;
2612 assert(compo < nbulles_tot);
2614 const double sf = surface_facette[fa7];
2617 for (
int dir = 0; dir < 3; dir++)
2621 for (
int j = 0; j < 3; j++)
2623 const int isom = facettes(fa7, j);
2624 v += vitesse_sommets(isom, dir);
2626 vitesses_translation_bulles(compo, dir) += v * sf / 3.;
2632 for (
int icompo = 0; icompo < nbulles_tot; icompo++)
2633 if (surface_par_bulle[icompo] > 0.)
2634 for (
int dir = 0; dir < 3; dir++)
2635 vitesses_translation_bulles(icompo, dir) /= surface_par_bulle[icompo];
2639 const ArrOfDouble& surface_facette,
2640 const ArrOfDouble& surface_par_bulle,
2641 const ArrOfInt& compo_connexes_facettes,
2642 const int nbulles_reelles,
2643 const int nbulles_ghost,
2644 const DoubleTab& centre_gravite,
2645 const DoubleTab& vitesse_sommets,
2646 const DoubleTab& vitesse_translation_sommets,
2647 DoubleTab& mean_bubble_rotation_vector)
const
2649 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2650 assert(mean_bubble_rotation_vector.
dimension(0) == nbulles_tot);
2651 assert(surface_par_bulle.
size_array() == nbulles_tot);
2652 assert(mean_bubble_rotation_vector.
dimension(1) == 3);
2654 const IntTab& facettes = maillage.
facettes();
2655 const DoubleTab& sommets = maillage.
sommets();
2657 mean_bubble_rotation_vector = 0.;
2664 for (
int fa7 = 0; fa7 < maillage.
nb_facettes(); fa7++)
2670 int compo = compo_connexes_facettes[fa7];
2678 compo = nbulles_reelles - 1 - idx_ghost;
2681 assert(compo < nbulles_tot);
2683 const double sf = surface_facette[fa7];
2686 for (
int j = 0; j < 3; j++)
2688 const int isom = facettes(fa7, j);
2690 double vit_x = vitesse_sommets(isom, 0) - vitesse_translation_sommets(compo, 0);
2691 double vit_y = vitesse_sommets(isom, 1) - vitesse_translation_sommets(compo, 1);
2692 double vit_z = vitesse_sommets(isom, 2) - vitesse_translation_sommets(compo, 2);
2693 Vecteur3 vit = {vit_x, vit_y, vit_z};
2695 double x = sommets(isom, 0);
2696 double y = sommets(isom, 1);
2697 double z = sommets(isom, 2);
2700 Vecteur3 coord_centre = {centre_gravite(compo,0), centre_gravite(compo,1), centre_gravite(compo,2)};
2702 Vecteur3 radial_position = coord_sommet - coord_centre;
2708 for (
int dir = 0; dir < 3; dir++)
2712 mean_bubble_rotation_vector(compo, dir) += rotation_vector[dir] * sf / 3.;
2719 for (
int icompo = 0; icompo < nbulles_tot; icompo++)
2720 if (surface_par_bulle[icompo] > 0.)
2721 for (
int dir = 0; dir < 3; dir++)
2722 mean_bubble_rotation_vector(icompo, dir) /= surface_par_bulle[icompo];
2734static void calculer_normale_sommets_interface(
const Maillage_FT_IJK& maillage,
2743 const IntTab& facettes = maillage.
facettes();
2744 const int nsommets_faces = facettes.
dimension(1);
2746 normale.
resize(nsom, dim);
2748 double n[3] = {0., 0. , 0.};
2750 ArrOfDouble surface_environnante(nsom);
2751 surface_environnante = 0.;
2753 for (
int i = 0; i < nfaces; i++)
2760 const double surface = surface_facettes[i];
2762 for (
int k = 0; k < dim; k++)
2763 n[k] = normale_facettes(i, k) * surface;
2765 for (
int j = 0; j < nsommets_faces; j++)
2767 const int sommet = facettes(i, j);
2768 surface_environnante[sommet] += surface;
2769 for (
int k = 0; k < dim; k++)
2770 normale(sommet, k) += n[k];
2781 int count[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2782 for (
int isom = 0; isom < nsom; isom++)
2785 for (
int dir = 0; dir < 3; dir++)
2787 if (surface_environnante[isom])
2788 normale(isom, dir) /= surface_environnante[isom];
2790 normale(isom, dir) = 0.;
2793 norme = normale(isom, 0) * normale(isom, 0) + normale(isom, 1) * normale(isom, 1) +
2794 normale(isom, 2) * normale(isom, 2);
2860 Cerr <<
"IJK_Interfaces.cpp:calculer_normale_sommets_interface : Calcul "
2861 "normale non-unitaire : ";
2863 for (
int i = 0; i < 9; i++)
2866 Cerr <<
" " << count[i] / nsom;
2868 Cerr <<
" " << 1 - sum / nsom <<
" (nsom = " << nsom << finl;
2873 const DoubleTab& vitesses_translation_bulles,
2874 const DoubleTab& mean_bubble_rotation_vector,
2875 const DoubleTab& centre_gravite,
2876 ArrOfDouble& var_volume)
2878 Cut_field_vector3_double& cut_field_deformation_velocity =
static_cast<Cut_field_vector3_double&
>(
deformation_velocity_);
2881 const DoubleTab& sommets = mesh.
sommets();
2886 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2889 var_volume.
resize(nbsom);
2901 for (
int dir = 0 ; dir < 3 ; dir++)
2906 DoubleTab bounding_box;
2912 for (
int icompo = 0; icompo < nbulles_tot; icompo++)
2922 cut_field_deformation_velocity,
2928 DoubleVect delta_volume_theorique_bilan_ns_vdf;
2932 const Domaine& domaine = domaine_vf.
domaine();
2938 assert(ni ==
I_ft().ni());
2939 assert(nj ==
I_ft().nj());
2940 assert(nk ==
I_ft().nk());
2942 for (
int k = 0; k < nk; k++)
2944 for (
int j = 0; j < nj; j++)
2946 for (
int i = 0; i < ni; i++)
2948 const int num_elem = ref_domaine_->convert_ijk_cell_to_packed(i, j, k);
2959 DoubleTab& mean_bubble_rotation_vector,
2960 DoubleTab& centre_gravite,
2961 const int first_step_interface_smoothing)
2967 if (first_step_interface_smoothing)
2980 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2988 ArrOfDouble surface_par_bulle;
2992 ArrOfDouble volume_par_bulle(nbulles_tot);
2993 vitesses_translation_bulles.
resize(nbulles_tot, 3);
2994 mean_bubble_rotation_vector.
resize(nbulles_tot, 3);
2995 centre_gravite.
resize(nbulles_tot, 3);
3010 vitesses_translation_bulles);
3025 vitesses_translation_bulles,
3026 mean_bubble_rotation_vector);
3028 if (first_step_interface_smoothing)
3030 vitesses_translation_bulles = 0.;
3031 mean_bubble_rotation_vector = 0.;
3040 for (
int bulle = 0; bulle < nbulles_tot; bulle++)
3042 Cerr <<
"composante " << bulle <<
" vitesse_translation ";
3043 for (
int i = 0; i < 3; i++)
3044 Cerr <<
" " << vitesses_translation_bulles(bulle, i);
3045 Cerr <<
" rotation_vector ";
3046 for (
int i = 0; i < 3; i++)
3047 Cerr <<
" " << mean_bubble_rotation_vector(bulle, i);
3065 const DoubleTab& vitesses_translation_bulles,
3066 const DoubleTab& mean_bubble_rotation_vector,
3067 const DoubleTab& centre_gravite,
3068 const double dt_tot,
3070 const int rk_step = -1,
3071 const int first_step_interface_smoothing)
3075 const DoubleTab& sommets = mesh.
sommets();
3077 DoubleTab deplacement(nbsom, 3);
3087 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
3095 ArrOfIntFT compo_connex_som;
3100 DoubleTabFT normale_sommet;
3101 calculer_normale_sommets_interface(mesh, normale_sommet);
3103 for (
int som = 0; som < nbsom; som++)
3108 deplacement(som, 0) = 100.;
3109 deplacement(som, 1) = 100.;
3110 deplacement(som, 2) = 100.;
3117 int icompo = compo_connex_som[som];
3125 icompo = nbulles_reelles - 1 - idx_ghost;
3127 assert(icompo >= 0);
3128 assert(icompo < nbulles_tot);
3146 double rot_x = mean_bubble_rotation_vector(icompo, 0);
3147 double rot_y = mean_bubble_rotation_vector(icompo, 1);
3148 double rot_z = mean_bubble_rotation_vector(icompo, 2);
3149 Vecteur3 rot = {rot_x, rot_y, rot_z};
3151 double x = sommets(som, 0);
3152 double y = sommets(som, 1);
3153 double z = sommets(som, 2);
3156 Vecteur3 coord_centre = {centre_gravite(icompo,0), centre_gravite(icompo,1), centre_gravite(icompo,2)};
3158 Vecteur3 radial_position = coord_sommet - coord_centre;
3163 double prodscal = 0.;
3164 double norme_carre = 0.;
3165 for (
int direction = 0; direction < 3; direction++)
3167 double vi =
vinterp_(som, direction);
3168 double vcompo = vitesses_translation_bulles(icompo, direction) + tangential_velocity[direction];
3169 double n = normale_sommet(som, direction);
3170 prodscal += (vi - vcompo) * n;
3171 norme_carre += n * n;
3174 prodscal /= norme_carre;
3178 for (
int direction = 0; direction < 3; direction++)
3180 double n = normale_sommet(som, direction);
3184 if (correction_semi_locale_volume_bulle)
3187 vinterp_(som, direction) = n * prodscal;
3192 vinterp_(som, direction) = n * prodscal + vitesses_translation_bulles(icompo, direction);
3221 Journal() <<
"RKSTEP " << rk_step << finl;
3222 for (
int i = 0; i < nbsom; i++)
3225 for (
int direction = 0; direction < 3; direction++)
3227 Journal() << vinterp(i, direction) <<
" ";
3229 Journal() <<
" RK3_G_store_vi_ ";
3230 for (
int direction = 0; direction < 3; direction++)
3252 for (
int som = 0; som < nbsom; som++)
3253 for (
int direction = 0; direction < 3; direction++)
3254 deplacement(som, direction) =
vinterp_(som, direction) * dt_tot;
3261 DoubleTab coord_sommets_avant_deplacement(mesh.
sommets());
3277 dt = compute_fractionnal_timestep_rk3(dt_tot, rk_step);
3321 if (!correction_semi_locale_volume_bulle)
3326 ArrOfDouble var_volume_par_bulle(nbulles_tot);
3331 for (
int isom = 0; isom < nbsom; isom++)
3337 int compo = compo_connex_som[isom];
3345 compo = nbulles_reelles - 1 - idx_ghost;
3348 assert(compo < nbulles_tot);
3351 var_volume_par_bulle[compo] += v;
3357 for (
int icompo = 0; icompo < nbulles_tot; icompo++)
3359 dvol[icompo] += var_volume_par_bulle[icompo];
3371 const DoubleTab& vitesses_translation_bulles,
3372 const DoubleTab& mean_bubble_rotation_vector,
3373 const DoubleTab& centre_gravite,
3376 const int rk_step = -1,
3377 const double temps = -1 )
3380 const DoubleTab& sommets = mesh.
sommets();
3383 if (correction_semi_locale_volume_bulle)
3385 const double fractionnal_timestep = (rk_step == -1) ? dt_tot : compute_fractionnal_timestep_rk3(dt_tot, rk_step);
3404 ArrOfDouble surface_par_bulle;
3420 if ((rk_step == -1) || (rk_step == 2))
3423 const IntTab& facettes = mesh.
facettes();
3424 for (
int i = 0; i < nb_facettes; i++)
3429 int compo = compo_connex_apres_transport[i];
3436 compo = nbulles_reelles - 1 - idx_ghost;
3438 const double var_volume_tot = dvol[compo];
3439 const double surface_bulle = surface_par_bulle[compo];
3440 const double s = surface_facette_apres_transport[i];
3441 const double dv = var_volume_tot * s / (surface_bulle * 3.);
3442 for (
int j = 0; j < 3; j++)
3444 const int isom = facettes(i, j);
3460 if ((rk_step == -1) || (rk_step == 2))
3486 const int nbsom_end_transport = mesh.
nb_sommets();
3488 if (!((nbsom_end_transport == size_store) || (0 == size_store)))
3490 Cerr <<
"Une des tailles de tableau n'est pas bonne... "
3491 <<
" size_store = " << size_store
3492 <<
" nbsom_end_transport = " << nbsom_end_transport
3506 const DoubleTab& vitesses_translation_bulles,
3507 const DoubleTab& mean_bubble_rotation_vector,
3508 const DoubleTab& centre_gravite,
3510 const int first_step_interface_smoothing)
3514 const DoubleTab& sommets = mesh.
sommets();
3516 DoubleTab deplacement(nbsom, 3);
3520 ArrOfIntFT compo_connex_som;
3523 for (
int som = 0; som < nbsom; som++)
3528 deplacement(som, 0) = 100.;
3529 deplacement(som, 1) = 100.;
3530 deplacement(som, 2) = 100.;
3534 int icompo = compo_connex_som[som];
3542 icompo = nbulles_reelles - 1 - idx_ghost;
3544 assert(icompo >= 0);
3547 double rot_x = mean_bubble_rotation_vector(icompo, 0);
3548 double rot_y = mean_bubble_rotation_vector(icompo, 1);
3549 double rot_z = mean_bubble_rotation_vector(icompo, 2);
3550 Vecteur3 rot = {rot_x, rot_y, rot_z};
3552 double x = sommets(som, 0);
3553 double y = sommets(som, 1);
3554 double z = sommets(som, 2);
3557 Vecteur3 coord_centre = {centre_gravite(icompo,0), centre_gravite(icompo,1), centre_gravite(icompo,2)};
3559 Vecteur3 radial_position = coord_sommet - coord_centre;
3564 for (
int direction = 0; direction < 3; direction++)
3566 const double fractionnal_timestep = (rk_step == -1) ? dt_tot : compute_fractionnal_timestep_rk3(dt_tot, rk_step);
3567 deplacement(som, direction) = (vitesses_translation_bulles(icompo, direction) + tangential_velocity[direction]) * fractionnal_timestep;
3577 ArrOfDouble& var_volume,
3596 if (algo_remaillage_local.
a_remailler(temps, maillage))
3608 const DoubleTab& bounding_box_bulles,
3609 const Cut_field_vector3_double& cut_field_velocity,
3610 const DoubleTab& vitesses_translation_bulles,
3611 const DoubleTab& mean_bubble_rotation_vector,
3612 const DoubleTab& positions_bulles)
3614 Cut_field_vector3_double& cut_field_deformation_velocity =
static_cast<Cut_field_vector3_double&
>(
deformation_velocity_);
3615 assert(&cut_field_deformation_velocity[0].get_cut_cell_disc() == &cut_field_deformation_velocity[1].get_cut_cell_disc());
3616 assert(&cut_field_deformation_velocity[0].get_cut_cell_disc() == &cut_field_deformation_velocity[2].get_cut_cell_disc());
3617 const Cut_cell_FT_Disc& cut_cell_disc = cut_field_deformation_velocity[0].get_cut_cell_disc();
3619 const int ni = cut_field_deformation_velocity[0].ni();
3620 const int nj = cut_field_deformation_velocity[0].nj();
3621 const int nk = cut_field_deformation_velocity[0].nk();
3622 const int ghost = cut_field_deformation_velocity[0].
ghost();
3633 double origin_x = geom.
get_origin(DIRECTION_I);
3634 double origin_y = geom.
get_origin(DIRECTION_J);
3635 double origin_z = geom.
get_origin(DIRECTION_K);
3641 int imin = std::max(-ghost, (
int)((bounding_box_bulles(compo, 0, 0) - origin_x)/dx - offset_x - 2));
3642 int imax = std::min(ni + ghost, (
int)((bounding_box_bulles(compo, 0, 1) - origin_x)/dx - offset_x + 2));
3643 int jmin = std::max(-ghost, (
int)((bounding_box_bulles(compo, 1, 0) - origin_y)/dy - offset_y - 2));
3644 int jmax = std::min(nj + ghost, (
int)((bounding_box_bulles(compo, 1, 1) - origin_y)/dy - offset_y + 2));
3645 int kmin = std::max(-ghost, (
int)((bounding_box_bulles(compo, 2, 0) - origin_z)/dz - offset_z - 2));
3646 int kmax = std::min(nk + ghost, (
int)((bounding_box_bulles(compo, 2, 1) - origin_z)/dz - offset_z + 2));
3647 for (
int k = kmin; k < kmax; k++)
3649 for (
int j = jmin; j < jmax; j++)
3651 for (
int i = imin; i < imax; i++)
3653 const double x_centre_cell = (i + offset_x + .5)*dx + origin_x;
3654 const double y_centre_cell = (j + offset_y + .5)*dy + origin_y;
3655 const double z_centre_cell = (k + offset_z + .5)*dz + origin_z;
3657 double rot_x = mean_bubble_rotation_vector(compo, 0);
3658 double rot_y = mean_bubble_rotation_vector(compo, 1);
3659 double rot_z = mean_bubble_rotation_vector(compo, 2);
3660 Vecteur3 rot = {rot_x, rot_y, rot_z};
3662 int n = cut_cell_disc.
get_n(i,j,k);
3665 for (
int phase = 0 ; phase < 2 ; phase++)
3671 Vecteur3 coord = {x_cut_cell, y_cut_cell, z_cut_cell};
3673 Vecteur3 coord_centre = {positions_bulles(compo,0), positions_bulles(compo,1), positions_bulles(compo,2)};
3675 Vecteur3 radial_position = coord - coord_centre;
3680 for (
int direction = 0; direction < 3; direction++)
3682 const DoubleTabFT_cut_cell& diph_velocity = (phase == 0) ? cut_field_velocity[direction].diph_v_ : cut_field_velocity[direction].diph_l_;
3683 DoubleTabFT_cut_cell& deformation_diph_velocity = (phase == 0) ? cut_field_deformation_velocity[direction].diph_v_ : cut_field_deformation_velocity[direction].diph_l_;
3684 deformation_diph_velocity(n) = diph_velocity(n) - vitesses_translation_bulles(compo, direction) - tangential_velocity[direction];
3690 Vecteur3 coord = {x_centre_cell, y_centre_cell, z_centre_cell};
3692 Vecteur3 coord_centre = {positions_bulles(compo,0), positions_bulles(compo,1), positions_bulles(compo,2)};
3694 Vecteur3 radial_position = coord - coord_centre;
3699 for (
int direction = 0; direction < 3; direction++)
3701 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];
3716 bounding_box.
resize(nbulles, 3, 2);
3720 const double inval = 1.e20;
3721 bounding_box = inval;
3725 ArrOfIntFT compo_connex_som;
3727 const DoubleTab& sommets = mesh.
sommets();
3729 ArrOfDouble volume_reel;
3733 ArrOfDouble position_xmax_compo;
3734 ArrOfDouble position_xmin_compo;
3737 position_xmax_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
3738 position_xmin_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
3739 position_xmax_compo = -10000.;
3740 position_xmin_compo = 10000.;
3742 for (
int i_sommet2 = 0; i_sommet2 < nbsom; i_sommet2++)
3744 int iconnex = compo_connex_som[i_sommet2];
3745 double coord = sommets(i_sommet2, 0);
3746 if (coord>position_xmax_compo[iconnex])
3747 position_xmax_compo[iconnex] = coord;
3748 if (coord<position_xmin_compo[iconnex])
3749 position_xmin_compo[iconnex] = coord;
3755 for (
int i_sommet = 0; i_sommet < nbsom; i_sommet++)
3757 for (direction=0; direction<3; direction++)
3759 double coord = sommets(i_sommet, direction);
3760 int iconnex = compo_connex_som[i_sommet];
3765 double pos_ref = position(iconnex,0);
3774 double decallage_bulle_reel_ext_domaine_reel = 1.*(position_xmax_compo[iconnex]-position_xmin_compo[iconnex]);
3777 double pos = std::fmod(std::fmod(pos_ref + offset - decallage_bulle_reel_ext_domaine_reel, Lx) + Lx, Lx) + decallage_bulle_reel_ext_domaine_reel;
3779 coord += (pos - pos_ref);
3786 double mmin = bounding_box(iconnex, direction, 0);
3787 double mmax = -bounding_box(iconnex, direction, 1);
3789 bounding_box(iconnex, direction, 0) = coord;
3791 bounding_box(iconnex, direction, 1) = -coord;
3802 for (
int ibulle = 0; ibulle < nbulles; ibulle++)
3803 for (direction = 0; direction < 3; direction++)
3804 bounding_box(ibulle, direction, 1) = -bounding_box(ibulle, direction, 1);
3813 DoubleTab bounding_box;
3815 ArrOfInt masque_duplicata_pour_compo_reel;
3828 DoubleTab bounding_box_offsetp;
3831 DoubleTab bounding_box_offsetm;
3850static int decoder_deplacement(
const int code,
const int dir,
int& compo_bulle_reel)
3853 const int code_deplacement = code & 63;
3854 int tmp = code_deplacement & (1 << (3 + dir));
3860 int signe = (tmp == 0) ? -1 : 1 ;
3861 tmp = code_deplacement & (1 << dir);
3863 int index = tmp >> dir;
3867 compo_bulle_reel = code >> 6;
3871 return signe * index;
3878 DoubleTab& deplacement,
3879 DoubleTab& bounding_box_NS, DoubleTab position,
3884 deplacement.
resize(nbsom, 3);
3886 ArrOfIntFT compo_sommets;
3890 ArrOfIntFT compo_sommets_init;
3894 ArrOfDouble position_xmax_compo;
3895 ArrOfDouble position_xmin_compo;
3898 position_xmax_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
3899 position_xmin_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
3900 position_xmax_compo = -10000.;
3901 position_xmin_compo = 10000.;
3903 const DoubleTab& sommets = mesh.
sommets();
3905 for (
int i_sommet = 0; i_sommet < nbsomreel; i_sommet++)
3909 int iconnex = compo_sommets_init[i_sommet];
3911 double coord = sommets(i_sommet, 0);
3912 if (coord>position_xmax_compo[iconnex])
3913 position_xmax_compo[iconnex] = coord;
3914 if (coord<position_xmin_compo[iconnex])
3915 position_xmin_compo[iconnex] = coord;
3921 for (
int dir = 0; dir < 3; dir++)
3923 for (
int i_sommet = 0; i_sommet < nbsom; i_sommet++)
3926 const int code = compo_sommets[i_sommet];
3931 int compo_bulle_reel = 0;
3933 int decode = decoder_deplacement(code, dir, compo_bulle_reel);
3935 double depl = decode * (bounding_box_NS(dir, 1) - bounding_box_NS(dir, 0));
3936 deplacement(i_sommet, dir) = depl;
3937 double pos_ref = 0 ;
3939 double decallage_bulle_reel_ext_domaine_reel = 0.;
3946 pos_ref = position(compo_bulle_reel,0);
3948 decallage_bulle_reel_ext_domaine_reel = 1.*(position_xmax_compo[compo_bulle_reel]-position_xmin_compo[compo_bulle_reel]);
3949 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;
3951 deplacement(i_sommet, 0) += (pos - pos_ref);
3963static void calculer_deplacement_from_code_compo_connexe_negatif(
const Maillage_FT_IJK& m,
3964 DoubleTab& deplacement,
3965 DoubleTab& bounding_box_NS)
3969 deplacement.
resize(nbsom, 3);
3970 ArrOfIntFT compo_sommets;
3973 for (
int i_sommet = 0; i_sommet < nbsom; i_sommet++)
3976 int code = compo_sommets[i_sommet];
3981 for (
int dir = 0; dir < 3; dir++)
3982 deplacement(i_sommet, dir) = 0.;
3990 for (
int dir = 0; dir < 3; dir++)
3993 int unused_variable = 0 ;
3994 int decode = decoder_deplacement(code, dir, unused_variable);
3995 double depl = decode * (bounding_box_NS(dir, 1) - bounding_box_NS(dir, 0));
3996 deplacement(i_sommet, dir) = depl;
4006static void calculer_deplacement_from_masque_in_array(
const Maillage_FT_IJK& m,
4007 DoubleTab& deplacement,
4008 const ArrOfInt& masque_array,
4009 DoubleTab& bounding_box_NS, DoubleTab position,
const int nbulles)
4013 deplacement.
resize(nbsom,3);
4014 ArrOfIntFT compo_connexe_som;
4017 ArrOfDouble position_xmax_compo;
4018 ArrOfDouble position_xmin_compo;
4021 position_xmax_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
4022 position_xmin_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
4023 position_xmax_compo = -10000.;
4024 position_xmin_compo = 10000.;
4026 const DoubleTab& sommets = m.
sommets();
4028 for (
int i_sommet = 0; i_sommet < nbsom; i_sommet++)
4030 int iconnex = compo_connexe_som[i_sommet];
4031 double coord = sommets(i_sommet, 0);
4032 if (coord>position_xmax_compo[iconnex])
4033 position_xmax_compo[iconnex] = coord;
4034 if (coord<position_xmin_compo[iconnex])
4035 position_xmin_compo[iconnex] = coord;
4041 for (
int i_sommet = 0; i_sommet < nbsom; i_sommet++)
4044 const int icompo = compo_connexe_som[i_sommet];
4045 const int code = masque_array[icompo];
4047 for (
int dir = 0; dir < 3; dir++)
4050 int compo_bulle_reel = 0;
4052 int decode = decoder_deplacement(code, dir, compo_bulle_reel);
4054 double depl = decode * (bounding_box_NS(dir, 1) - bounding_box_NS(dir, 0));
4055 deplacement(i_sommet, dir) = depl;
4056 double pos_ref = 0 ;
4058 double decallage_bulle_reel_ext_domaine_reel = 0.;
4065 pos_ref = position(compo_bulle_reel,0);
4067 decallage_bulle_reel_ext_domaine_reel = 1.*(position_xmax_compo[compo_bulle_reel]-position_xmin_compo[compo_bulle_reel]);
4068 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;
4070 deplacement(i_sommet, 0) += (pos - pos_ref);
4095 maillage_temporaire.
initialize(ref_domaine_.valeur(), refdomaine_dis_.valeur(),
parcours_);
4106 int icompo, index, signe;
4107 for (
int i_facette = 0; i_facette < a_dupliquer.
nb_facettes(); i_facette++)
4109 icompo = compo_connexe_facettes[i_facette];
4126 ArrOfInt liste_bulles_crees;
4128 int nbulles_crees = 0;
4129 for (
int mon_numero_iteration = 1; ; mon_numero_iteration++)
4133 ArrOfInt index_copie(nbulles);
4134 ArrOfInt index_signe(nbulles);
4137 ArrOfInt liste_facettes_pour_suppression;
4138 int reste_a_faire = 0;
4142 for (icompo = 0; icompo < nbulles; icompo++)
4146 index_copie[icompo] = -1;
4150 const int masque = masque_duplicata_pour_compo[icompo] & 7;
4151 for (index = 1; index <= 7; index++)
4153 if ((masque & index) == index)
4156 if (compteur == mon_numero_iteration)
4160 index_copie[icompo] = index;
4171 index_signe[icompo] = -1;
4172 ArrOfInt index_bulle(7);
4173 ArrOfInt signe_bulle(7);
4175 const int perio_x_reel = masque_duplicata_pour_compo[icompo] & 1;
4176 const int perio_y_reel = masque_duplicata_pour_compo[icompo] & (1 << 1);
4177 const int perio_z_reel = masque_duplicata_pour_compo[icompo] & (1 << 2);
4178 const int perio_x_ghost = masque_duplicata_pour_compo[icompo] & (1 << 3);
4181 int signe_3bit = (masque_duplicata_pour_compo[icompo] & (7 << 4)) >> 1;
4186 int bit_position_x_signe_6bit = 3;
4187 int bit_position_x_ghost_8bit = 7;
4188 int signe_x_ghost = (masque_duplicata_pour_compo[icompo] & (1 << bit_position_x_ghost_8bit)) >> bit_position_x_ghost_8bit;
4189 int mask = 1<<bit_position_x_signe_6bit;
4191 for (
int nb_bulle_duplique = 0; nb_bulle_duplique < 7; nb_bulle_duplique++)
4193 index_bulle[nb_bulle_duplique] = -1;
4194 signe_bulle[nb_bulle_duplique] = signe_3bit;
4216 if (perio_z_reel == 4)
4218 if (perio_y_reel == 2 && perio_x_reel == 1 && perio_x_ghost == 8 )
4226 signe_bulle[4] &= (~mask);
4227 signe_bulle[4] |= (signe_x_ghost << bit_position_x_signe_6bit);
4230 signe_bulle[6] &= (~mask);
4231 signe_bulle[6] |= (signe_x_ghost << bit_position_x_signe_6bit);
4233 else if (perio_y_reel == 2 && perio_x_reel == 1 && perio_x_ghost != 8 )
4241 else if (perio_y_reel == 2 && perio_x_reel != 1 && perio_x_ghost == 8 )
4244 signe_bulle[0] &= (~mask);
4245 signe_bulle[0] |= (signe_x_ghost << bit_position_x_signe_6bit);
4248 signe_bulle[2] &= (~mask);
4249 signe_bulle[2] |= (signe_x_ghost << bit_position_x_signe_6bit);
4253 else if (perio_y_reel != 2 && perio_x_reel == 1 && perio_x_ghost == 8 )
4258 signe_bulle[2] &= (~mask);
4259 signe_bulle[2] |= (signe_x_ghost << bit_position_x_signe_6bit);
4262 else if (perio_y_reel != 2 && perio_x_reel != 1 && perio_x_ghost != 8 )
4266 else if (perio_y_reel != 2 && perio_x_reel != 1 && perio_x_ghost == 8 )
4269 signe_bulle[0] &= (~mask);
4270 signe_bulle[0] |= (signe_x_ghost << bit_position_x_signe_6bit);
4273 else if (perio_y_reel != 2 && perio_x_reel == 1 && perio_x_ghost != 8 )
4278 else if (perio_y_reel == 2 && perio_x_reel != 1 && perio_x_ghost != 8 )
4286 Cerr <<
"Erreur dans la gestion des bulles ghost" << finl;
4292 if (perio_y_reel == 2 && perio_x_reel == 1)
4298 else if (perio_y_reel == 2 && perio_x_reel != 1)
4302 else if (perio_y_reel != 2 && perio_x_reel == 1)
4306 else if (perio_y_reel != 2 && perio_x_reel != 1)
4308 index_bulle[0] = -1;
4312 Cerr <<
"Erreur dans la gestion des bulles ghost" << finl;
4317 if (index_bulle[mon_numero_iteration-1] != -1 && mon_numero_iteration-1 <= 6)
4320 index_copie[icompo] = index_bulle[mon_numero_iteration-1];
4321 index_signe[icompo] = signe_bulle[mon_numero_iteration-1];
4328 if (reste_a_faire == 0)
4336 for (
int i_facette = 0; i_facette < nf; i_facette++)
4338 icompo = decoder_numero_bulle(compo_connexe_facettes[i_facette]);
4340 index = index_copie[icompo];
4342 signe = index_signe[icompo];
4365 signe = masque_duplicata_pour_compo[icompo] & (7 << 3);
4366 const int code_deplacement = signe | index;
4377 Cerr <<
"Comment est-ce possible? " << finl;
4382 liste_facettes_pour_suppression.
append_array(i_facette);
4406 DoubleTab deplacement;
4411 ArrOfDouble volume_reel;
4415 calculer_deplacement_from_code_compo_connexe(maillage_temporaire, split,
4427 for (
int iface = 0; iface < nbf; iface++)
4429 const int code_negatif = -compo_connexe_facettes_tempo[iface];
4439 array_trier_retirer_doublons(liste_bulles_crees);
4450 for (
int p = 1; p < nbproc; p++)
4452 recevoir(prov, p, 0, 2001 + p);
4462 array_trier_retirer_doublons(liste_bulles_crees);
4465 nbulles_crees = liste_bulles_crees.
size_array();
4472 Cerr <<
" On vient de creer des bulles ghost alors qu'il y en avait "
4474 <<
" tableau ghost_compo_converter_ n'est pas de taille nulle! " << finl;
4491 Cerr <<
"IJK_Interfaces.cpp::dupliquer_bulle_perio liste_bulles_crees : "
4492 << liste_bulles_crees << finl;
4496 if (new_size != old_size)
4498 Cerr <<
"On dirait que des doublons etaient presents dans la liste ghost_compo_converter_ ... "
4499 <<
"Ca signifie peut-etre que l'on vient de creer des bulles qui existaient deja..."
4500 <<
"C'est etrange... on s'arrete... ";
4534 const DoubleTab& bounding_box_offsetp,
4535 const DoubleTab& bounding_box_offsetm,
4536 const DoubleTab& authorized_bounding_box,
4537 ArrOfInt& masque_duplicata_pour_compo)
4544 for (
int icompo = 0; icompo < nbulles; icompo++)
4546 int masque_sortie_domaine_reel = 0;
4552 for (
int direction = 0; direction < 3; direction++)
4557 if (bounding_box(icompo, direction, 0) < authorized_bounding_box(direction, 0))
4559 masque_sortie_domaine_reel |= (1 << direction);
4560 masque_sortie_domaine_reel |= (16 << direction);
4567 if (bounding_box_offsetp(icompo, 0, 0) < authorized_bounding_box(0, 0))
4569 masque_sortie_domaine_reel |= (1 << 3);
4570 masque_sortie_domaine_reel |= (16 << 3);
4572 if (bounding_box_offsetp(icompo, 0, 1) > authorized_bounding_box(0, 1))
4574 masque_sortie_domaine_reel |= (1 << 3);
4580 if (bounding_box(icompo, direction, 1) > authorized_bounding_box(direction, 1))
4584 masque_sortie_domaine_reel |= (1 << direction);
4592 if (bounding_box_offsetm(icompo, 0, 0) < authorized_bounding_box(0, 0))
4594 masque_sortie_domaine_reel |= (1 << 3);
4595 masque_sortie_domaine_reel |= (16 << 3);
4597 if (bounding_box_offsetm(icompo, 0, 1) > authorized_bounding_box(0, 1))
4599 masque_sortie_domaine_reel |= (1 << 3);
4609 masque_duplicata_pour_compo[icompo] = masque_sortie_domaine_reel;
4614 envoyer_broadcast(masque_duplicata_pour_compo, 0);
4619 const DoubleTab& authorized_bounding_box,
4620 ArrOfInt& masque_duplicata_pour_compo)
4627 for (
int icompo = 0; icompo < nbulles; icompo++)
4629 int masque_sortie_domaine_reel = 0;
4635 for (
int direction = 0; direction < 3; direction++)
4640 if (bounding_box(icompo, direction, 0) < authorized_bounding_box(direction, 0))
4642 masque_sortie_domaine_reel |= (1 << direction);
4643 masque_sortie_domaine_reel |= (8 << direction);
4646 if (bounding_box(icompo, direction, 1) > authorized_bounding_box(direction, 1))
4648 masque_sortie_domaine_reel |= (1 << direction);
4654 masque_duplicata_pour_compo[icompo] = masque_sortie_domaine_reel;
4660 envoyer_broadcast(masque_duplicata_pour_compo, 0);
4673 ArrOfInt liste_facettes_pour_suppression;
4674 for (
int iface = 0; iface < nbfacettes; iface++)
4676 const int icompo = compo_connex[iface];
4696 DoubleTab bounding_box;
4714 for (
int icompo = 0; icompo < nbulles; icompo++)
4716 const int masque_sortie = masque[icompo];
4717 if (masque_sortie & 2)
4721 if (masque_sortie & 16)
4752 DoubleTab deplacement;
4754 ArrOfDouble volume_reel;
4757 calculer_deplacement_from_masque_in_array(mesh, deplacement, masque_deplacement_par_compo,
bounding_box_NS_domain_, position, nbulles);
4763 for (
int ib = 0; ib < nbulles; ib++)
4765 const int code_deplacement = masque_deplacement_par_compo[ib];
4766 if (code_deplacement != 0)
4768 Journal() <<
"IJK_Interfaces::deplacer_bulle_perio : Un deplacement a eu "
4769 <<
"lieu pour la composante " << ib << finl;
4770 Journal() <<
"Deplacement x y z : ";
4771 for (
int dir = 0; dir < 3; dir++)
4773 int unused_variable = 0 ;
4774 const int decode = decoder_deplacement(code_deplacement, dir, unused_variable);
4807 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul de l'indicatrice",2,
"IJK");
4808 statistics().begin_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice",statistics().get_last_opened_counter_level()+1);
4811 const IntTab& elem_faces = domaine_vf.
elem_faces();
4812 const IntTab& faces_voisins = domaine_vf.
face_voisins();
4818 const int ni = indic.
ni();
4819 const int nj = indic.
nj();
4820 const int nk = indic.
nk();
4828 const ArrOfInt& index_elem = intersec.
index_elem();
4831 for (
int k = 0; k < nk; k++)
4833 for (
int j = 0; j < nj; j++)
4835 for (
int i = 0; i < ni; i++)
4842 int index = index_elem[num_elem];
4843 double somme_contrib = 0.;
4852 int nb_increment_somme_contrib = 0;
4853 while (somme_contrib > 1.)
4855 somme_contrib -= 1.;
4856 nb_increment_somme_contrib++;
4858 while (somme_contrib < 0.)
4860 somme_contrib += 1.;
4861 nb_increment_somme_contrib++;
4864 if (nb_increment_somme_contrib > 1)
4866 Cerr <<
"nb_increment_somme_contrib= " << nb_increment_somme_contrib << finl;
4878 if (somme_contrib > 0.)
4883 indic(i, j, k) = somme_contrib;
4891 statistics().create_custom_counter(
"Calcul de l'indicatrice : recherche compo connexes",2,
"IJK");
4892 statistics().begin_count(
"Calcul de l'indicatrice : recherche compo connexes",statistics().get_last_opened_counter_level()+1);
4895 int nb_compo_locales = search_connex_components_local(elem_faces, faces_voisins,
num_compo_);
4902 statistics().end_count(
"Calcul de l'indicatrice : recherche compo connexes");
4909 Cerr <<
"GB-check-interf nb_compo_in_num_compo_ - (nb_reelles + nb_ghost)= " <<
nb_compo_in_num_compo_ <<
" - "
4913 Cerr <<
"IJK_Interfaces::calculer_indicatrice Evaluations drapeau vapeur " << drapeau << finl;
4916#define VERIF_DRAPEAU 0
4931 compute_drapeaux_vapeur_v1(mesh, elem_faces, faces_voisins, num_compo,
true, drapeau1);
4932 compute_drapeaux_vapeur_v1(mesh, elem_faces, faces_voisins, num_compo,
false, drapeau1b);
4941 Cerr <<
"IJK_Interfaces::calculer_indicatrice Evaluations drapeaux" << finl;
4943 Cerr <<
"drapeau v1 towards_liquid=true : " << drapeau1 << finl;
4944 Cerr <<
"drapeau v1 towards_liquid=false : " << drapeau1b << finl;
4945 Cerr <<
"drapeau v2 : " << drapeau2 << finl;
4946 Cerr <<
"drapeau v3 : " << drapeau3 << finl;
4947 Cerr <<
"drapeau v4 : " << drapeau << finl;
4951 int sum = drapeau[idx] + drapeau1[idx] + drapeau1b[idx] + drapeau2[idx] + drapeau3[idx];
4952 if ((drapeau[idx] != drapeau1[idx]) || (drapeau[idx] != drapeau1b[idx]) || (drapeau[idx] != drapeau2[idx]) ||
4953 (drapeau[idx] != drapeau3[idx]))
4954 Cerr <<
"Les methodes different... "
4955 <<
"valeurs donnees par les methodes : " << drapeau[idx] << drapeau1[idx] << drapeau1b[idx]
4956 << drapeau2[idx] << drapeau3[idx] <<
" pour idx = " << idx << finl;
4958 if ((sum != 0) && (sum != 5))
4960 Cerr <<
"Les methodes different... "
4961 <<
"sum = " << sum <<
" pour idx = " << idx << finl;
4965 assert(drapeau == drapeau1);
4966 assert(drapeau == drapeau1b);
4967 assert(drapeau == drapeau2);
4968 assert(drapeau == drapeau3);
4971 check_somme_drapeau(drapeau);
4972 check_somme_drapeau(drapeau1);
4973 check_somme_drapeau(drapeau1b);
4974 check_somme_drapeau(drapeau2);
4975 check_somme_drapeau(drapeau3);
4979#define DEBUG_INDIC 0
4984 for (
int k = 0; k < nk; k++)
4986 for (
int j = 0; j < nj; j++)
4988 for (
int i = 0; i < ni; i++)
4991 int compo = num_compo[num_elem];
4993 indic(i, j, k) = compo;
5002 for (
int k = 0; k < nk; k++)
5004 for (
int j = 0; j < nj; j++)
5006 for (
int i = 0; i < ni; i++)
5011 if (compo >= 0 && drapeau[compo] == 0)
5014 indic(i, j, k) = 0.;
5021 statistics().end_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice");
5026 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul de l'indicatrice",2,
"IJK");
5027 statistics().begin_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice",statistics().get_last_opened_counter_level()+1);
5032 const int ni = indic.
ni();
5033 const int nj = indic.
nj();
5034 const int nk = indic.
nk();
5039 const ArrOfInt& index_elem = intersec.
index_elem();
5041 for (
int k = 0; k < nk; k++)
5043 for (
int j = 0; j < nj; j++)
5045 for (
int i = 0; i < ni; i++)
5048 int index = index_elem[num_elem];
5050 if (0. < indic(i, j, k) && indic(i, j, k) < 1.)
5051 indic(i, j, k) = 2.;
5053 double somme_contrib = 0.;
5062 int nb_increment_somme_contrib = 0;
5063 while (somme_contrib > 1.)
5065 somme_contrib -= 1.;
5066 nb_increment_somme_contrib++;
5068 while (somme_contrib < 0.)
5070 somme_contrib += 1.;
5071 nb_increment_somme_contrib++;
5074 if (nb_increment_somme_contrib > 1)
5075 Process::exit(
"Error in IJK_Interfaces::calculer_indicatrice_optim !");
5077 if (somme_contrib > 0.)
5081 indic(i, j, k) = somme_contrib;
5089 statistics().end_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice");
5094 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul des indicatrices",2,
"IJK");
5095 statistics().begin_count(
"Calcul rho mu indicatrice: calcul des indicatrices",statistics().get_last_opened_counter_level()+1);
5097 const Domaine& domaine = domaine_vf.
domaine();
5098 const IntTab& elem_faces = domaine_vf.
elem_faces();
5099 const IntTab& faces_voisins = domaine_vf.
face_voisins();
5103 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5105 indic[igroup].data() = 1.;
5109 Cerr <<
"This choice of options has to be validated. "
5110 <<
"It seems you are using a multi-groups of bubbles "
5111 <<
"and the forcing method with the color_function "
5112 <<
"which has not been tested yet!"
5117 IntVect num_compo[max_authorized_nb_of_groups_];
5122 domaine.creer_tableau_elements(num_compo[i]);
5124 const int ni = indic[0].ni();
5125 const int nj = indic[0].nj();
5126 const int nk = indic[0].nk();
5128 assert(
nb_groups_ <= max_authorized_nb_of_groups_);
5137 const ArrOfInt& index_elem = intersec.
index_elem();
5140 for (
int k = 0; k < nk; k++)
5142 for (
int j = 0; j < nj; j++)
5144 for (
int i = 0; i < ni; i++)
5151 int index = index_elem[num_elem];
5152 double somme_contrib[max_authorized_nb_of_groups_] = {0.};
5158 int icompo = compo_connexe[facette];
5164 m = decoder_numero_bulle(-icompo);
5172 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5174 while (somme_contrib[igroup] > 1.)
5175 somme_contrib[igroup] -= 1.;
5176 while (somme_contrib[igroup] < 0.)
5177 somme_contrib[igroup] += 1.;
5178 if (somme_contrib[igroup] > 0.)
5182 indic[igroup](i, j, k) = somme_contrib[igroup];
5183 num_compo[igroup][num_elem] = -1;
5191 ArrOfInt drapeau[max_authorized_nb_of_groups_];
5192 statistics().create_custom_counter(
"Calcul de l'indicatrice : recherche compo connexes",2,
"IJK");
5195 statistics().begin_count(
"Calcul de l'indicatrice : recherche compo connexes",statistics().get_last_opened_counter_level()+1);
5198 int nb_compo_locales = search_connex_components_local(elem_faces, faces_voisins, num_compo[i]);
5201 statistics().end_count(
"Calcul de l'indicatrice : recherche compo connexes");
5211 for (
int k = 0; k < nk; k++)
5213 for (
int j = 0; j < nj; j++)
5215 for (
int i = 0; i < ni; i++)
5218 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5220 int compo = num_compo[igroup][num_elem];
5221 if (compo >= 0 && drapeau[igroup][compo] == 0)
5224 indic[igroup](i, j, k) = 0.;
5232 statistics().end_count(
"Calcul rho mu indicatrice: calcul des indicatrices");
5237 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul des indicatrices",2,
"IJK");
5238 statistics().begin_count(
"Calcul rho mu indicatrice: calcul des indicatrices",statistics().get_last_opened_counter_level()+1);
5243 const int ni = indic[0].ni();
5244 const int nj = indic[0].nj();
5245 const int nk = indic[0].nk();
5247 assert(
nb_groups_ <= max_authorized_nb_of_groups_);
5253 const ArrOfInt& index_elem = intersec.
index_elem();
5256 for (
int k = 0; k < nk; k++)
5258 for (
int j = 0; j < nj; j++)
5260 for (
int i = 0; i < ni; i++)
5263 int index = index_elem[num_elem];
5264 double somme_contrib[max_authorized_nb_of_groups_] = {0.};
5266 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5268 if (0. < indic[igroup](i, j, k) && indic[igroup](i, j, k) < 1.)
5269 indic[igroup](i, j, k) = 2.;
5277 int icompo = compo_connexe[facette];
5283 m = decoder_numero_bulle(-icompo);
5290 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5292 while (somme_contrib[igroup] > 1.)
5293 somme_contrib[igroup] -= 1.;
5294 while (somme_contrib[igroup] < 0.)
5295 somme_contrib[igroup] += 1.;
5296 if (somme_contrib[igroup] > 0.)
5298 indic[igroup](i, j, k) = somme_contrib[igroup];
5306 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5310 statistics().end_count(
"Calcul rho mu indicatrice: calcul des indicatrices");
5315 const int ni = indic.
ni();
5316 const int nj = indic.
nj();
5317 const int nk = indic.
nk();
5337 for (
int k = 0; k < nk; k++)
5339 for (
int j = 0; j < nj; j++)
5341 for (
int i = 0; i < ni; i++)
5343 if (indic(i, j, k) == 2.)
5345 ArrOfInt voisins_num(6);
5353 ArrOfDouble voisins_valeur(6);
5354 voisins_valeur[0] = imin + i == 0 ? BORD : indic(i - 1, j, k);
5355 voisins_valeur[1] = jmin + j == 0 ? BORD : indic(i, j - 1, k);
5356 voisins_valeur[2] = kmin + k == 0 ? BORD : indic(i, j, k - 1);
5357 voisins_valeur[3] = imin + i == nitot ? BORD : indic(i + 1, j, k);
5358 voisins_valeur[4] = jmin + j == njtot ? BORD : indic(i, j + 1, k);
5359 voisins_valeur[5] = kmin + k == nktot ? BORD : indic(i, j, k + 1);
5366 for (
int v = 0; v < 6; v++)
5368 double valeur_indic_voisin = voisins_valeur[v];
5369 if (valeur_indic_voisin == 0. || valeur_indic_voisin == 1.)
5371 indic(i, j, k) = valeur_indic_voisin;
5383 for (
int v = 0; v < 6; v++)
5385 double valeur_indic_voisin = voisins_valeur[v];
5386 const int direction = v % 3;
5387 const int face_plus = (v > 2) ? 1 : -1;
5389 if (0. < valeur_indic_voisin && valeur_indic_voisin < 1.)
5391 int num_elem = voisins_num[v];
5392 if (num_elem != BORD)
5397 if (indic(i, j, k) != -1)
5429 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face",2,
"IJK");
5430 statistics().begin_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face",statistics().get_last_opened_counter_level()+1);
5435 const int ni = indic_surfacique_face[0].ni();
5436 const int nj = indic_surfacique_face[0].nj();
5437 const int nk = indic_surfacique_face[0].nk();
5441 for (
int k = 0; k < nk; k++)
5443 for (
int j = 0; j < nj; j++)
5445 for (
int i = 0; i < ni; i++)
5450 indic_surfacique_face[0](i, j, k) = indic(i, j, k);
5451 indic_surfacique_face[1](i, j, k) = indic(i, j, k);
5452 indic_surfacique_face[2](i, j, k) = indic(i, j, k);
5458 indic_surfacique_face[0](i, j, k) = norme[0](i, j, k) == 0. ? 0. : (norme[0](i, j, k) > 0. ? 0. : 1.);
5459 indic_surfacique_face[1](i, j, k) = norme[1](i, j, k) == 0. ? 0. : (norme[1](i, j, k) > 0. ? 0. : 1.);
5460 indic_surfacique_face[2](i, j, k) = norme[2](i, j, k) == 0. ? 0. : (norme[2](i, j, k) > 0. ? 0. : 1.);
5466 indic_surfacique_face[0](i, j, k) = indic(i-1, j, k);
5470 indic_surfacique_face[1](i, j, k) = indic(i, j-1, k);
5474 indic_surfacique_face[2](i, j, k) = indic(i, j, k-1);
5479 baric_face[0][0](i, j, k) = 0.5;
5480 baric_face[0][1](i, j, k) = 0.5;
5481 baric_face[1][0](i, j, k) = 0.5;
5482 baric_face[1][1](i, j, k) = 0.5;
5483 baric_face[2][0](i, j, k) = 0.5;
5484 baric_face[2][1](i, j, k) = 0.5;
5493 const ArrOfInt& index_elem = intersec.
index_elem();
5496 for (
int k = 0; k < nk; k++)
5498 for (
int j = 0; j < nj; j++)
5500 for (
int i = 0; i < ni; i++)
5512 int index = index_elem[num_elem];
5513 double somme_contrib[3] = {0., 0., 0.};
5514 double somme_contrib_baryc[3][2] = {{0., 0.}, {0., 0.}, {0., 0.}};
5533 for (
int dir=0; dir<3; dir++)
5535 const int dir_i = (dir == 0);
5536 const int dir_j = (dir == 1);
5537 const int dir_k = (dir == 2);
5539 if (
est_pure(indic(i-dir_i,j-dir_j,k-dir_k)))
5549 if (somme_contrib[dir]*somme_contrib[dir] > 0)
5551 while (somme_contrib[dir] > 1.)
5553 somme_contrib[dir] -= 1.;
5554 somme_contrib_baryc[dir][0] -= 1./2.;
5555 somme_contrib_baryc[dir][1] -= 1./2.;
5557 while (somme_contrib[dir] < 0.)
5559 somme_contrib[dir] += 1.;
5560 somme_contrib_baryc[dir][0] += 1./2.;
5561 somme_contrib_baryc[dir][1] += 1./2.;
5564 indic_surfacique_face[dir](i, j, k) = somme_contrib[dir];
5568 baric_face[dir][0](i, j, k) = 1./2.;
5569 baric_face[dir][1](i, j, k) = 1./2.;
5573 baric_face[dir][0](i, j, k) = somme_contrib_baryc[dir][0]/abs(somme_contrib[dir]);
5574 baric_face[dir][1](i, j, k) = somme_contrib_baryc[dir][1]/abs(somme_contrib[dir]);
5578 if ((baric_face[dir][0](i, j, k) <= 0) && (baric_face[dir][0](i, j, k) >= -1e-12 || somme_contrib[dir] < 1e-8))
5580 baric_face[dir][0](i, j, k) = 1e-12;
5582 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))
5584 baric_face[dir][0](i, j, k) = 1 - 1e-12;
5586 if ((baric_face[dir][1](i, j, k) <= 0) && (baric_face[dir][1](i, j, k) >= -1e-12 || somme_contrib[dir] < 1e-8))
5588 baric_face[dir][1](i, j, k) = 1e-12;
5590 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))
5592 baric_face[dir][1](i, j, k) = 1 - 1e-12;
5595 assert((baric_face[dir][0](i, j, k) >= 0) && (baric_face[dir][0](i, j, k) <= 1));
5596 assert((baric_face[dir][1](i, j, k) >= 0) && (baric_face[dir][1](i, j, k) <= 1));
5597 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.)));
5598 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.)));
5608 statistics().end_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face");
5614 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face",2,
"IJK");
5615 statistics().begin_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face",statistics().get_last_opened_counter_level()+1);
5620 const int ni = indic_surfacique_face[0].ni();
5621 const int nj = indic_surfacique_face[0].nj();
5622 const int nk = indic_surfacique_face[0].nk();
5626 for (
int k = 0; k < nk; k++)
5628 for (
int j = 0; j < nj; j++)
5630 for (
int i = 0; i < ni; i++)
5635 indic_surfacique_face[0](i, j, k) = indic(i, j, k);
5636 indic_surfacique_face[1](i, j, k) = indic(i, j, k);
5637 indic_surfacique_face[2](i, j, k) = indic(i, j, k);
5643 indic_surfacique_face[0](i, j, k) = norme[0](i, j, k) == 0. ? 0. : (norme[0](i, j, k) > 0. ? 0. : 1.);
5644 indic_surfacique_face[1](i, j, k) = norme[1](i, j, k) == 0. ? 0. : (norme[1](i, j, k) > 0. ? 0. : 1.);
5645 indic_surfacique_face[2](i, j, k) = norme[2](i, j, k) == 0. ? 0. : (norme[2](i, j, k) > 0. ? 0. : 1.);
5651 indic_surfacique_face[0](i, j, k) = indic(i-1, j, k);
5655 indic_surfacique_face[1](i, j, k) = indic(i, j-1, k);
5659 indic_surfacique_face[2](i, j, k) = indic(i, j, k-1);
5670 const ArrOfInt& index_elem = intersec.
index_elem();
5673 for (
int k = 0; k < nk; k++)
5675 for (
int j = 0; j < nj; j++)
5677 for (
int i = 0; i < ni; i++)
5689 int index = index_elem[num_elem];
5690 double somme_contrib[3] = {0., 0., 0.};
5702 for (
int dir=0; dir<3; dir++)
5704 const int dir_i = (dir == 0);
5705 const int dir_j = (dir == 1);
5706 const int dir_k = (dir == 2);
5708 if (
est_pure(indic(i-dir_i,j-dir_j,k-dir_k)))
5718 if (somme_contrib[dir]*somme_contrib[dir] > 0)
5720 while (somme_contrib[dir] > 1.)
5722 somme_contrib[dir] -= 1.;
5724 while (somme_contrib[dir] < 0.)
5726 somme_contrib[dir] += 1.;
5729 indic_surfacique_face[dir](i, j, k) = somme_contrib[dir];
5738 statistics().end_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face");
5744 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul de la surface interfaciale",2,
"IJK");
5745 statistics().begin_count(
"Calcul rho mu indicatrice: calcul de la surface interfaciale",statistics().get_last_opened_counter_level()+1);
5749 const ArrOfDouble& surface_facettes =
maillage_ft_ijk_.get_update_surface_facettes();
5751 const int ni = surf_interface.
ni();
5752 const int nj = surf_interface.
nj();
5753 const int nk = surf_interface.
nk();
5757 for (
int k = 0; k < nk; k++)
5759 for (
int j = 0; j < nj; j++)
5761 for (
int i = 0; i < ni; i++)
5763 surf_interface(i, j, k) = 0.;
5772 const ArrOfInt& index_elem = intersec.
index_elem();
5775 for (
int k = 0; k < nk; k++)
5777 for (
int j = 0; j < nj; j++)
5779 for (
int i = 0; i < ni; i++)
5790 int index = index_elem[num_elem];
5791 double somme_contrib_surf = 0.;
5802 surf_interface(i, j, k) = somme_contrib_surf;
5803 assert(somme_contrib_surf >= 0.);
5809 statistics().end_count(
"Calcul rho mu indicatrice: calcul de la surface interfaciale");
5814 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul du barycentre",2,
"IJK");
5815 statistics().begin_count(
"Calcul rho mu indicatrice: calcul du barycentre",statistics().get_last_opened_counter_level()+1);
5819 const int ni = baric[0].ni();
5820 const int nj = baric[0].nj();
5821 const int nk = baric[0].nk();
5825 for (
int k = 0; k < nk; k++)
5827 for (
int j = 0; j < nj; j++)
5829 for (
int i = 0; i < ni; i++)
5832 baric[0](i, j, k) = 0.5;
5833 baric[1](i, j, k) = 0.5;
5834 baric[2](i, j, k) = 0.5;
5843 const ArrOfInt& index_elem = intersec.
index_elem();
5846 for (
int k = 0; k < nk; k++)
5848 for (
int j = 0; j < nj; j++)
5850 for (
int i = 0; i < ni; i++)
5861 int index = index_elem[num_elem];
5862 double somme_contrib = 0.;
5863 double somme_contrib_baryc[3] = {0., 0., 0.};
5878 while (somme_contrib > 1.)
5880 somme_contrib -= 1.;
5881 somme_contrib_baryc[0] -= 1./2.;
5882 somme_contrib_baryc[1] -= 1./2.;
5883 somme_contrib_baryc[2] -= 1./2.;
5886 while (somme_contrib < 0.)
5888 somme_contrib += 1.;
5889 somme_contrib_baryc[0] += 1./2.;
5890 somme_contrib_baryc[1] += 1./2.;
5891 somme_contrib_baryc[2] += 1./2.;
5896 assert(indic(i, j, k) == somme_contrib);
5898 baric[0](i, j, k) = somme_contrib_baryc[0]/abs(somme_contrib);
5899 baric[1](i, j, k) = somme_contrib_baryc[1]/abs(somme_contrib);
5900 baric[2](i, j, k) = somme_contrib_baryc[2]/abs(somme_contrib);
5904 if ((baric[0](i, j, k) <= 0) && (baric[0](i, j, k) >= -1e-12 || somme_contrib < 1e-8))
5906 baric[0](i, j, k) = 1e-12;
5908 else if ((baric[0](i, j, k) >= 1) && (baric[0](i, j, k) <= 1+1e-12 || somme_contrib < 1e-8))
5910 baric[0](i, j, k) = 1 - 1e-12;
5912 if ((baric[1](i, j, k) <= 0) && (baric[1](i, j, k) >= -1e-12 || somme_contrib < 1e-8))
5914 baric[1](i, j, k) = 1e-12;
5916 else if ((baric[1](i, j, k) >= 1) && (baric[1](i, j, k) <= 1+1e-12 || somme_contrib < 1e-8))
5918 baric[1](i, j, k) = 1 - 1e-12;
5920 if ((baric[2](i, j, k) <= 0) && (baric[2](i, j, k) >= -1e-12 || somme_contrib < 1e-8))
5922 baric[2](i, j, k) = 1e-12;
5924 else if ((baric[2](i, j, k) >= 1) && (baric[2](i, j, k) <= 1+1e-12 || somme_contrib < 1e-8))
5926 baric[2](i, j, k) = 1 - 1e-12;
5929 assert((baric[0](i, j, k) >= 0) && (baric[0](i, j, k) <= 1));
5930 assert((baric[1](i, j, k) >= 0) && (baric[1](i, j, k) <= 1));
5931 assert((baric[2](i, j, k) >= 0) && (baric[2](i, j, k) <= 1));
5932 assert(((indic(i, j, k) < 0.) && (baric[0](i, j, k) < 0.)) || ((indic(i, j, k) > 0.) && (baric[0](i, j, k) > 0.)));
5933 assert(((indic(i, j, k) < 0.) && (baric[1](i, j, k) < 0.)) || ((indic(i, j, k) > 0.) && (baric[1](i, j, k) > 0.)));
5934 assert(((indic(i, j, k) < 0.) && (baric[2](i, j, k) < 0.)) || ((indic(i, j, k) > 0.) && (baric[2](i, j, k) > 0.)));
5940 statistics().end_count(
"Calcul rho mu indicatrice: calcul du barycentre");
5951 const Domaine& domaine = domaine_vf.
domaine();
5954 const int nb_elem = domaine.nb_elem();
5955 const int nb_elem_tot = domaine.nb_elem_tot();
5956 Cerr <<
"nb_elem= " << nb_elem <<
" nb_elem_tot= " << nb_elem_tot;
5971#define NUM_COMPO_INVALID (-2000000000)
5979 IJK_Field_vector3_double& vpoint,
5980 IJK_Field_vector3_double& vrepul,
5981 IJK_Field_vector3_double& vabsrepul
5984 statistics().begin_count(STD_COUNTERS::source_terms,statistics().get_last_opened_counter_level()+1);
5996 double interf1, interf2 ;
5998 Int3 ijk_other_elem;
5999 for (
int icol1 = 0; icol1 < max_authorized_nb_of_components_; icol1++)
6001 for (
int direction = 0; direction < 3; direction++)
6003 for (
int k = 0; k < vpoint[direction].nk(); k++)
6005 for (
int j = 0; j < vpoint[direction].nj(); j++)
6007 for (
int i = 0; i < vpoint[direction].ni(); i++)
6014 ijk_other_elem[0] = i;
6015 ijk_other_elem[1] = j;
6016 ijk_other_elem[2] = k;
6017 ijk_other_elem[direction] -= 1;
6019 s2 =
surf_par_compo_[
old()][icol1](ijk_other_elem[0],ijk_other_elem[1],ijk_other_elem[2]);
6022 vpoint[direction](i, j, k) += 1./2. * (s1 * interf1 + s2 * interf2);
6035 const int nkmax = std::max(vpoint[DIRECTION_I].nk(), std::max(vpoint[DIRECTION_J].nk(), vpoint[DIRECTION_K].nk()));
6036 for (
int k = 0; k < nkmax; k++)
6038 for (
int direction = 0; direction < 3; direction++)
6040 if (k >= vpoint[direction].nk())
6048 for (
int j = 0; j < vpoint[direction].nj(); j++)
6050 for (
int i = 0; i < vpoint[direction].ni(); i++)
6055 const int global_face_position = ijk_face[direction] + offset;
6056 if (!perio && (global_face_position == 0 || global_face_position == nb_items_tot))
6059 Int3 ijk_droite = ijk_face;
6060 Int3 ijk_gauche = ijk_face;
6061 ijk_gauche[direction]--;
6065 for (
int gauche_droite = 0; gauche_droite <= 1; gauche_droite++)
6068 for (
int icol1 = 0; icol1 < max_authorized_nb_of_components_; icol1++)
6070 const Int3& elem1 = gauche_droite ? ijk_droite : ijk_gauche;
6071 const Int3& elem2 = gauche_droite ? ijk_gauche : ijk_droite;
6073 int signe = gauche_droite ? -1 : 1;
6075 if (num_compo == NUM_COMPO_INVALID)
6080 for (icol2 = 0; icol2 < max_authorized_nb_of_components_; icol2++)
6084 if (icol2 < max_authorized_nb_of_components_ && gauche_droite == 1)
6097 elem1[0], elem1[1], elem1[2]
6103 double indic_voisin = 0., phi_voisin = 0.;
6104 double repul_voisin = 0.;
6105 if (icol2 == max_authorized_nb_of_components_)
6136 for (
int dir = 0; dir < 3; dir++)
6138 const int idx = icol1 * 3 + dir;
6145 bary_facettes_dans_elem[dir] =
bary_par_compo_[
old()][idx](elem1[0], elem1[1], elem1[2]);
6152 indic_voisin = (ps > 0 ? 1. : 0.);
6155 repul_voisin = repul;
6169 double gradient_phi_indic = (phi_voisin * indic_voisin - phi * indic) / delta_dir * signe;
6171 vpoint[direction](i, j, k) += gradient_phi_indic;
6172 double gradient_repul_indic = (repul_voisin * indic_voisin - repul * indic) / delta_dir * signe;
6173 vrepul[direction](i, j, k) += gradient_repul_indic;
6174 vabsrepul[direction](i, j, k) += fabs(gradient_repul_indic);
6178 double indic_voisin = 0., phi_face = 0.;
6179 double repul_face = 0.;
6180 if (icol2 == max_authorized_nb_of_components_)
6211 for (
int dir = 0; dir < 3; dir++)
6213 const int idx = icol1 * 3 + dir;
6216 bary_facettes_dans_elem[dir] =
bary_par_compo_[
old()][idx](elem1[0], elem1[1], elem1[2]);
6223 indic_voisin = (ps > 0 ? 1. : 0.);
6241 elem2[0], elem2[1], elem2[2]
6246 assert(indic_voisin > -0.5);
6247 if (est_egal(surface+surface_voisin,0.))
6255 phi_face = (phi * surface + phi_voisin * surface_voisin) / (surface + surface_voisin);
6256 repul_face = (repul * surface + repul_voisin * surface_voisin) / (surface + surface_voisin);
6260 double gradient_indic = (indic_voisin - indic) / delta_dir * signe;
6270 vpoint[direction](i, j, k) += phi_face * gradient_indic ;
6271 vrepul[direction](i, j, k) += repul_face * gradient_indic;
6272 vabsrepul[direction](i, j, k) += fabs(repul_face * gradient_indic);
6281 statistics().end_count(STD_COUNTERS::source_terms);
6293static bool intersection_segment_triangle(
const Vecteur3& A,
6307 const double det1 = determinant(CD, CE, CA);
6308 const double det2 = determinant(CD, CE, CB);
6310 if (det1 * det2 >= 0)
6320 const double det3 = determinant(AC, AD, AB);
6321 const double det4 = determinant(AD, AE, AB);
6322 const double det5 = determinant(AE, AC, AB);
6324 if ((det3 * det4 <= 0) || (det3 * det5 <= 0))
6336 const IntTab& facettes = mesh.
facettes();
6337 const DoubleTab& sommets = mesh.
sommets();
6342 ArrOfInt facettes_traversantes;
6344 const int N = facettes_traversantes.
size_array();
6353 int fa7 = facettes_traversantes[index];
6373 for (
int i = 0; i < 3; i++)
6374 sommets_facette[i] =
Vecteur3(sommets, facettes(fa7, i));
6382 for (
int isom = 0; isom < 3; isom++)
6385 for (
int dir = 0; dir < 3; dir++)
6386 coordA[dir] += bary_som * sommets_facette[isom][dir];
6393 coordB[direction] += cell_size[direction] * face_plus;
6398 for (
int index2 = intersections.
index_elem()[num_elem]; index2 >= 0;
6406 const Vecteur3 s0(sommets, facettes(nvl_fa7, 0));
6407 const Vecteur3 s1(sommets, facettes(nvl_fa7, 1));
6408 const Vecteur3 s2(sommets, facettes(nvl_fa7, 2));
6411 coupe = intersection_segment_triangle(coordA, coordB, s0, s1, s2);
6417 for (
int j = 0; j < N; j++)
6419 if (facettes_traversantes[j] == nvl_fa7)
6434 if (normale_facette[direction] * face_plus > 0.)
6447#define COMPUTE_DRAPEAUX_VALIDATION
6450 ArrOfInt& drapeau_vapeur)
const
6452 statistics().create_custom_counter(
"Calcul de l'indicatrice: calculs des drapeaux",2,
"IJK");
6453 statistics().begin_count(
"Calcul de l'indicatrice: calculs des drapeaux",statistics().get_last_opened_counter_level()+1);
6458 const ArrOfInt& index_elem =
maillage_ft_ijk_.intersections_elem_facettes().index_elem();
6459 const ArrOfInt& index_facette =
maillage_ft_ijk_.intersections_elem_facettes().index_facette();
6462 const IntTab& elem_faces = domaine_vf.
elem_faces();
6463 const IntTab& faces_voisins = domaine_vf.
face_voisins();
6465#ifdef COMPUTE_DRAPEAUX_VALIDATION
6466 ArrOfInt composantes_comptage_phase0(drapeau_vapeur.
size_array());
6467 ArrOfInt composantes_comptage_phase1(drapeau_vapeur.
size_array());
6476 const IntTab& facettes = mesh.
facettes();
6477 const DoubleTab& sommets = mesh.
sommets();
6478 for (
int fa7 = 0; fa7 < nb_facettes; fa7++)
6485 for (
int i = 0; i < 3; i++)
6486 sommets_facette[i] =
Vecteur3(sommets, facettes(fa7, i));
6488 int index = index_facette[fa7];
6516 for (
int iface = 0; iface < 6; iface++)
6518 const int direction = iface % 3;
6519 const int face_plus = (iface > 2) ? 1 : -1;
6521 const int num_face = elem_faces(elem, iface);
6522 const int elem_voisin = faces_voisins(num_face, 0) + faces_voisins(num_face, 1) - elem;
6524 if (elem_voisin < 0)
6530 const int global_compo_voisin = vecteur_composantes[elem_voisin];
6534 if (global_compo_voisin < 0)
6537#ifndef COMPUTE_DRAPEAUX_VALIDATION
6540 if (drapeau_vapeur[global_compo_voisin] == 1)
6545 sommets_facette[2] - sommets_facette[0],
6548#ifndef COMPUTE_DRAPEAUX_VALIDATION
6552 if (normale_facette[direction] * face_plus <= 0.)
6561 for (
int isom = 0; isom < 3; isom++)
6564 for (
int dir = 0; dir < 3; dir++)
6565 coordA[dir] += bary_som * sommets_facette[isom][dir];
6574 coordB[direction] += cell_size[direction] * face_plus;
6583 for (
int index2 = index_elem[elem]; index2 >= 0;
6591 const Vecteur3 s0(sommets, facettes(nvl_fa7, 0));
6592 const Vecteur3 s1(sommets, facettes(nvl_fa7, 1));
6593 const Vecteur3 s2(sommets, facettes(nvl_fa7, 2));
6596 coupe = intersection_segment_triangle(coordA, coordB, s0, s1, s2);
6607 if (normale_facette[direction] * face_plus > 0.)
6608 drapeau_vapeur[global_compo_voisin] = 1;
6609#ifdef COMPUTE_DRAPEAUX_VALIDATION
6611 if (normale_facette[direction] * face_plus > 0.)
6613 composantes_comptage_phase1[global_compo_voisin]++;
6621 if (normale_facette[direction] * face_plus < 0.)
6623 composantes_comptage_phase0[global_compo_voisin]++;
6643#ifdef COMPUTE_DRAPEAUX_VALIDATION
6645 const int n = composantes_comptage_phase0.
size_array();
6651 for (
int i = 0; i < n; i++)
6656 if ((composantes_comptage_phase0[i] > 0 && composantes_comptage_phase1[i] > 0) ||
6657 (composantes_comptage_phase0[i] == 0 && composantes_comptage_phase1[i] == 0))
6662 Cerr <<
"[WARNING-Indic-Vote] compo/#phase0/#phase1 ";
6665 <<
"/" << composantes_comptage_phase0[i]
6666 <<
"/" << composantes_comptage_phase1[i]
6670 if (composantes_comptage_phase0[i] > composantes_comptage_phase1[i])
6671 drapeau_vapeur[i] = 0;
6673 drapeau_vapeur[i] = 1;
6677 Cerr <<
"End." << finl;
6679 envoyer_broadcast(drapeau_vapeur, 0);
6683 statistics().end_count(
"Calcul de l'indicatrice: calculs des drapeaux");
6686static inline double norme_carre(
const Vecteur3& x)
6688 return x[0] * x[0] + x[1] * x[1] + x[2] * x[2];
6693 return Vecteur3(x * v[0], x * v[1], x * v[2]);
6697 return Vecteur3(v[0] + w[0], v[1] + w[1], v[2] + w[2]);
6700static inline double calculer_carre_distance_sommet_facette(
const Vecteur3& coord,
6704 const IntTab& facettes = mesh.
facettes();
6705 const DoubleTab& sommets = mesh.
sommets();
6706 Vecteur3 s0(sommets, facettes(num_facette, 0));
6707 Vecteur3 s1(sommets, facettes(num_facette, 1));
6708 Vecteur3 s2(sommets, facettes(num_facette, 2));
6715 for (
double x = 0.; x < 1.01; x += 0.2)
6716 for (
double y = 0.; y < 1. - x + 0.01; y += 0.2)
6718 Vecteur3 s = s0 + x * (s1 - s0) + y * (s2 - s0);
6719 d = std::min(d, norme_carre(s - coord));
6724static void fill_relative_velocity(
const DoubleTab& vinterp_tmp,
const DoubleTab& vinterp,
const IntTab& facettes,
int id_facette,
int som, DoubleTab& vr_to_other)
6726 const double un_tiers = 1. / 3.;
6727 if (id_facette == -1)
6731 for (
int idir = 0; idir < 3; idir++)
6732 vr_to_other(som, idir) = 0.;
6737 const int isom0 = facettes(id_facette, 0);
6738 const int isom1 = facettes(id_facette, 1);
6739 const int isom2 = facettes(id_facette, 2);
6740 for (
int idir = 0; idir < 3; idir++)
6744 const double velocity_me = vinterp_tmp(som, idir);
6745 const double velocity_other =
6746 un_tiers * (vinterp(isom0, idir) + vinterp(isom1, idir) + vinterp(isom2, idir));
6747 vr_to_other(som, idir) = velocity_me - velocity_other;
6755 const ArrOfInt& compo_connexe_sommets,
6756 const DoubleTab& vinterp_tmp,
6758 ArrOfDouble& distance,
6759 DoubleTab& vr_to_other,
6760 const double distmax)
6763 Octree_Double octree;
6764 const IntTab& facettes = mesh.
facettes();
6769 ArrOfInt liste_facettes;
6772 const int nb_som = sommets_a_tester.
dimension(0);
6773 distance.resize_array(nb_som, RESIZE_OPTIONS::NOCOPY_NOINIT);
6775 vr_to_other.
resize(nb_som, 3, RESIZE_OPTIONS::NOCOPY_NOINIT);
6777 assert(vinterp_tmp.
dimension(0) == nb_som);
6778 for (
int i = 0; i < nb_som; i++)
6780 Vecteur3 coord(sommets_a_tester, i);
6781 const int compo_connexe_som = compo_connexe_sommets[i];
6783 double dmin = distmax * distmax;
6786 octree.
search_elements_box(coord[0] - distmax, coord[1] - distmax, coord[2] - distmax, coord[0] + distmax,
6787 coord[1] + distmax, coord[2] + distmax, liste_facettes);
6791 const int nliste = liste_facettes.
size_array();
6792 int idx_facette = -1;
6793 for (
int j = 0; j < nliste; j++)
6795 const int num_facette = liste_facettes[j];
6796 const int compo = compo_connexe_facettes[num_facette];
6797 if (compo == compo_connexe_som)
6801 const double d = calculer_carre_distance_sommet_facette(coord, mesh, num_facette);
6805 idx_facette = num_facette;
6809 distance[i] = std::min(distance[i], sqrt(dmin));
6810 fill_relative_velocity(vinterp_tmp,
vinterp_, facettes, idx_facette, i, vr_to_other);
6817static void check_neighbouring_layer_in_one_direction(
int dir0,
int dir1,
int dir2,
int n,
6818 const Int3& nb_elem_loc,
const Int3& ijk,
6819 const std::map<std::array<int,3>, std::set<int>>& bary_ijk_loc,
const DoubleTab& bary,
6820 const ArrOfInt& compo_connexe_facettes,
const int compo_connexe_som,
6821 const double x,
const double y,
const double z,
6822 double& distance,
int& id_facette )
6825 for(
int sens=0; sens<2; sens++)
6828 int a = sens == 0 ? std::max(ijk[dir0]-n,zero) : std::min(ijk[dir0]+n,nb_elem_loc[dir0]);
6829 for(
int b=std::max(ijk[dir1]-n,zero); b<std::min(ijk[dir1]+n,nb_elem_loc[dir1]); b++)
6830 for(
int c=std::max(ijk[dir2]-n,zero); c<std::min(ijk[dir2]+n,nb_elem_loc[dir2]); c++)
6832 std::array<int,3> current_ijk;
6836 if(bary_ijk_loc.find(current_ijk)!=bary_ijk_loc.end())
6838 std::set<int> bary_in_current_ijk = bary_ijk_loc.at(current_ijk);
6839 for(
const auto b_fa7: bary_in_current_ijk)
6841 if(compo_connexe_facettes[b_fa7] != compo_connexe_som)
6844 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);
6845 distance = std::min(sqrt(dist),distance);
6873 const ArrOfInt& compo_connexe_sommets,
6874 const DoubleTab& vinterp_tmp,
6876 ArrOfDouble& distance,
6877 DoubleTab& vr_to_other,
6878 const double distmax)
6880 const IntTab& facettes = mesh.
facettes();
6882 const int nb_som = sommets_a_tester.
dimension(0);
6883 const int nb_fa7 = facettes.
dimension(0);
6886 Int3 ijk_glob, ijk_loc, useless;
6888 for(
int dir=0; dir<3; dir++)
6891 distance.resize_array(nb_som, RESIZE_OPTIONS::NOCOPY_NOINIT);
6893 vr_to_other.
resize(nb_som, 3, RESIZE_OPTIONS::NOCOPY_NOINIT);
6895 assert(vinterp_tmp.
dimension(0) == nb_som);
6897 DoubleTab bary(nb_fa7,3);
6899 std::map<std::array<int,3>, std::set<int>> bary_ijk_loc;
6900 for (
int fa7=0; fa7<nb_fa7; fa7++)
6903 int s0 = facettes(fa7,0), s1 = facettes(fa7,1), s2 = facettes(fa7,2);
6904 for(
int dir=0; dir<3; dir++)
6905 bary(fa7,dir) = (sommets_a_tester(s0,dir) + sommets_a_tester(s1,dir) + sommets_a_tester(s2,dir)) / 3.;
6907 splitting.
search_elem(bary(fa7,0), bary(fa7,1), bary(fa7,2), ijk_glob, ijk_loc, useless);
6908 std::array<int,3> ijk;
6909 for(
int dir=0; dir<3; dir++) ijk[dir] = ijk_loc[dir];
6910 bary_ijk_loc[ijk].insert(fa7);
6913 for (
int som=0; som<nb_som; som++)
6915 double x=sommets_a_tester(som,0), y=sommets_a_tester(som,1), z=sommets_a_tester(som,2);
6916 splitting.
search_elem(x, y, z, ijk_glob, ijk_loc, useless);
6917 const int compo_connexe_som = compo_connexe_sommets[som];
6919 bool local_domain_checked =
false;
6921 int id_facette = -1;
6922 double& dist = distance[som];
6925 while(!local_domain_checked && id_facette==-1)
6928 check_neighbouring_layer_in_one_direction(0 , 1, 2, n_layer,
6929 nb_elem_loc, ijk_loc,
6931 compo_connexe_facettes, compo_connexe_som,
6936 check_neighbouring_layer_in_one_direction(1 , 0, 2, n_layer,
6937 nb_elem_loc, ijk_loc,
6939 compo_connexe_facettes, compo_connexe_som,
6946 check_neighbouring_layer_in_one_direction(2 , 0, 1, n_layer,
6947 nb_elem_loc, ijk_loc,
6949 compo_connexe_facettes, compo_connexe_som,
6955 int i=ijk_loc[0], j=ijk_loc[1], k=ijk_loc[2];
6956 local_domain_checked = i-n_layer <= 0 && i+n_layer>= nb_elem_loc[0]
6957 && j-n_layer <= 0 && j+n_layer>= nb_elem_loc[1]
6958 && k-n_layer <= 0 && k+n_layer>= nb_elem_loc[2];
6965 fill_relative_velocity(vinterp_tmp,
vinterp_, facettes, id_facette, som, vr_to_other);
6992 ArrOfInt& compo_sommet, ArrOfDouble& distance,
6993 DoubleTab& vr_to_other,
double distmax)
6995 const Domaine_IJK& splitting = ref_domaine_.valeur();
7010 double min_coord, max_coord;
7012 if (processor_at_left < 0)
7020 const double left_node = coord_nodes[offset];
7021 min_coord = left_node + distmax;
7022 if(!flags[processor_at_left])
7024 flags[processor_at_left] = 1;
7029 if (processor_at_right < 0)
7038 const double right_node = coord_nodes[i_last_node];
7039 max_coord = right_node - distmax;
7040 if(!flags[processor_at_right])
7042 flags[processor_at_right] = 1;
7049 ArrOfInt index_sent_to_left;
7050 ArrOfInt index_sent_to_right;
7051 if (processor_at_left >= 0 || processor_at_right >= 0)
7054 const int nb_som = coord_sommets.
dimension(0);
7055 for (
int i_som = 0; i_som < nb_som; i_som++)
7057 Vecteur3 coord(coord_sommets, i_som);
7058 if (coord[dir] < min_coord)
7064 schema.
send_buffer(processor_at_left) << compo_sommet[i_som]
7065 << coord_sommets(i_som, 0)
7066 << coord_sommets(i_som, 1)
7067 << coord_sommets(i_som, 2)
7068 << vinterp_tmp(i_som, 0)
7069 << vinterp_tmp(i_som, 1)
7070 << vinterp_tmp(i_som, 2);
7072 if (coord[dir] > max_coord)
7076 schema.
send_buffer(processor_at_right) << compo_sommet[i_som]
7077 << coord_sommets(i_som, 0)
7078 << coord_sommets(i_som, 1)
7079 << coord_sommets(i_som, 2)
7080 << vinterp_tmp(i_som, 0)
7081 << vinterp_tmp(i_som, 1)
7082 << vinterp_tmp(i_som, 2);
7089 const int init_count = coord_sommets.
dimension(0);
7091 int count_right = 0;
7092 for (
int left_right = 0; left_right < 2; left_right++)
7094 const int pe_voisin = (left_right==0) ? processor_at_left : processor_at_right;
7095 int& count = (left_right==0) ? count_left : count_right;
7105 double x, y, z, vx, vy, vz;
7106 buf >> x >> y >> z >> vx >> vy >> vz;
7124 int index = init_count;
7126 for (
int left_right = 0; left_right < 2; left_right++)
7128 const int pe_voisin = (left_right==0) ? processor_at_left : processor_at_right;
7129 const int count = (left_right==0) ? count_left : count_right;
7133 for (
int i = 0; i < count; i++)
7134 buf << distance[index+i]
7135 << vr_to_other(index+i,0)
7136 << vr_to_other(index+i,1)
7137 << vr_to_other(index+i,2);
7142 for (
int left_right = 0; left_right < 2; left_right++)
7144 int pe_voisin = (left_right==0) ? processor_at_left : processor_at_right;
7145 const ArrOfInt indices_sommets = (left_right==0) ? index_sent_to_left : index_sent_to_right;
7149 const int count = indices_sommets.
size_array();
7150 for (
int i = 0; i < count; i++)
7155 buf >> vx >> vy >> vz;
7156 int idx = indices_sommets[i];
7157 distance[idx] = std::min(distance[idx], d);
7158 vr_to_other(idx,0) = vx;
7159 vr_to_other(idx,1) = vy;
7160 vr_to_other(idx,2) = vz;
7167 coord_sommets.
resize(init_count, 3);
7168 vinterp_tmp.
resize(init_count, 3);
7171 distance.resize_array(init_count);
7172 vr_to_other.
resize(init_count, 3);
7181 DoubleTab& v_closer)
7185 ArrOfIntFT compo_connexe_sommets;
7187 DoubleTab tmp_sommets = mesh.
sommets();
7197 compo_connexe_sommets,
7210 f.
ouvrir(
"dump_facette.txt", ios::app);
7215 f.
ouvrir(
"dump_facette.txt");
7217 for (
int i = 0; i < 4; i++)
7219 int som = mesh.
facettes()(fa7, i%3);
7220 f << mesh.
sommets()(som, 0) <<
" " << mesh.
sommets()(som, 1) <<
" " << mesh.
sommets()(som, 2) << finl;
7224void dump_elem(
const Domaine_IJK& split,
int ii,
int jj,
int kk,
int append = 0)
7229 f.
ouvrir(
"dump_elem.txt", ios::app);
7234 f.
ouvrir(
"dump_elem.txt");
7241 for (
int i = 0; i < 3; i++)
7247 for (
int dir = 0; dir < 3; dir++)
7249 for (
int seg = 0; seg < 8; seg++)
7252 int pend = seg + (1 << dir);
7255 f << p[pstart & 1][0] <<
" " << p[(pstart >> 1) & 1][1] <<
" " << p[(pstart >> 2) & 1][2] << finl;
7256 f << p[pend & 1][0] <<
" " << p[(pend >> 1) & 1][1] <<
" " << p[(pend >> 2) & 1][2] << finl;
7263void dump_elem(
const Domaine_VF& domaine,
int elem,
int append = 0)
7268 f.
ouvrir(
"dump_elem.txt", ios::app);
7273 f.
ouvrir(
"dump_elem.txt");
7276 for (
int i = 0; i < 3; i++)
7278 p[0][i] = domaine.xv(domaine.elem_faces(elem, i), i);
7279 p[1][i] = domaine.xv(domaine.elem_faces(elem, i+3), i);
7281 for (
int dir = 0; dir < 3; dir++)
7283 for (
int seg = 0; seg < 8; seg++)
7286 int pend = seg + (1 << dir);
7289 f << p[pstart & 1][0] <<
" " << p[(pstart >> 1) & 1][1] <<
" " << p[(pstart >> 2) & 1][2] << finl;
7290 f << p[pend & 1][0] <<
" " << p[(pend >> 1) & 1][1] <<
" " << p[(pend >> 2) & 1][2] << finl;
7311 if (duplicatas_etaient_presents)
7323 ArrOfInt compo_new(nb_facettes);
7324 int n = search_connex_components_local_FT(maillage, compo_new);
7326 int nb_compo_tot=compute_global_connex_components_FT(maillage, compo_new, n);
7331 ArrOfInt rejeton(nb_rejeton);
7338 for (
int fa7 = 0; fa7 < nb_facettes; fa7++)
7339 rejetonArea(compo_std[fa7], compo_new[fa7])+=surface_facettes[fa7];
7349 for (
int ibul_new = 0 ; ibul_new<nb_compo_tot ; ibul_new++)
7350 if (rejetonArea(ibul, ibul_new)!=0.0 && rejetonArea(ibul, ibul_stay)<rejetonArea(ibul, ibul_new))
7352 ibul_stay=ibul_new ;
7357 for (
int ibul_new = 0 ; ibul_new<nb_compo_tot ; ibul_new++)
7358 if (rejetonArea(ibul, ibul_new)!=0.0 && ibul_new!=ibul_stay)
7359 for (
int ii = 0; ii<nb_rejeton ; ii++)
7360 if (rejeton[ii]==-1)
7362 rejeton[ii]=ibul_new;
7366 Cerr <<
"matrice de surface des compo" << finl;
7367 Cerr << rejetonArea << finl;
7368 Cerr <<
"Les rejetons a supprimer sont" << finl;
7369 Cerr << rejeton << finl;
7372 for (
int fa7 = 0; fa7 < nb_facettes; fa7++)
7373 for (
int ii = 0; ii<nb_rejeton ; ii++)
7374 if (compo_new[fa7] == rejeton[ii])
7382 if (duplicatas_etaient_presents)
7396 IJK_Field_vector3_double& rappel,
7397 const IJK_Field_vector3_double& vitesse,
7398 const IJK_Field_double& indic,
7399 const IJK_Field_double& indic_ft,
7400 const double coef_immo,
7402 const double current_time,
7403 const double coef_ammortissement,
7404 const double coef_rayon_force_rappel,
7405 const double integration_time,
7406 const double coef_mean_force,
7407 const double coef_force_time_n)
7426 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
7430 ArrOfDouble surface_par_bulle;
7433 ArrOfIntFT compo_connex_som;
7436 DoubleTab vitesses_translation_bulles(nbulles_tot,3);
7444 vitesses_translation_bulles);
7449 ArrOfDouble volume_reel;
7453 const double deltat = 1.;
7455 for (
int idir=0; idir < 3; idir++)
7464 while (dx> 0.8*ldom)
7467 while (dx< -0.8*ldom)
7470 individual_forces(ib,idir) = coef_mean_force*
mean_force_(ib,idir);
7471 individual_forces(ib,idir) += (1.-coef_mean_force)*coef_immo*(dx);
7472 individual_forces(ib,idir) += (1.-coef_mean_force)*coef_ammortissement*vitesses_translation_bulles(ib,idir);
7473 individual_forces(ib,idir) += (1.-coef_mean_force)*coef_force_time_n*
force_time_n_(ib,idir);
7485 for (
int idir=0; idir < 3; idir++)
7489 mean_force_(ib,idir)+=individual_forces(ib,idir)*deltat;
7494 envoyer_broadcast(individual_forces, 0);
7520 int reset = (!
reprise_) && (tstep==0);
7521 IOS_OPEN_MODE mode = (reset) ? ios::out : ios::app;
7522 for (
int idir=0; idir < 3; idir++)
7524 snprintf(s, 1000,
"%s_bulles_external_force_every_%d.out", nomcas, idir);
7530 snprintf(s, 1000,
"# Individual forces applied inside chiv[bubble_i]=1.\n# value=1./Vol_bubble \\int F_j(bubble_i) dv");
7534 snprintf(s, 1000,
"%.16e ", current_time);
7538 snprintf(s, 1000,
"%.16e ", individual_forces(ib,idir));
7551 const IJK_Field_double& indic_ns,
7552 const IJK_Field_double& indic_ft,
7553 DoubleTab& individual_forces,
7554 const ArrOfDouble& volume_reel,
7555 const DoubleTab& position)
7560 decodeur_num_compo=NUM_COMPO_INVALID;
7566 Int3 ijk_global, ijk_local, ijk_processeur;
7570 double x = position(ib,0);
7571 double y = position(ib,1);
7572 double z = position(ib,2);
7573 ijk_processeur[0]=0;
7574 ijk_processeur[1]=0;
7575 ijk_processeur[2]=0;
7578 s_ft.
search_elem(x,y,z, ijk_global, ijk_local, ijk_processeur);
7580 if (ijk_processeur[0] * ijk_processeur[1] * ijk_processeur[2] == 1)
7585 const int icompo =
num_compo_[num_elem_zvdf];
7587 decodeur_num_compo[icompo] = ib;
7595 ArrOfInt list_continuous_phase;
7597 int phase_continue1=-1000;
7598 int phase_continue2=-1000;
7599 int phase_continue3=-1000;
7600 int phase_continue4=-1000;
7601 int phase_continue5=-1000;
7602 int phase_continue6=-1000;
7606 if (decodeur_num_compo[i] == NUM_COMPO_INVALID)
7609 if (phase_continue1<0)
7611 phase_continue1 = i;
7613 else if (phase_continue2<0)
7615 phase_continue2 = i;
7617 else if (phase_continue3<0)
7619 phase_continue3 = i;
7621 else if (phase_continue4<0)
7623 phase_continue4 = i;
7625 else if (phase_continue5<0)
7627 phase_continue5 = i;
7629 else if (phase_continue6<0)
7631 phase_continue6 = i;
7636 const int nb_parts_continuous = list_continuous_phase.
size_array();
7637 Cerr <<
"nb parts continuous phase : " << list_continuous_phase.
size_array()
7638 <<
" continuous phases are : " << phase_continue1 <<
" " << phase_continue2 <<
" "
7639 << phase_continue3 <<
" " << phase_continue4 <<
" " << phase_continue5 <<
" " << phase_continue6 << finl;
7640 Cerr <<
"GB-check nb_compo - nb_continuous phase - nb_bulles_tot = "
7645 if (nb_parts_continuous == 0)
7647 Cerr <<
"No continuous phase found. " << finl;
7663 integration_cells_per_bubble=0;
7664 integrated_forces =0.;
7665 for (
int idir=0; idir < 3; idir++)
7667 const int nx = rappel_ft[idir].ni();
7668 const int ny = rappel_ft[idir].nj();
7670 int nz = rappel_ft[idir].nk();
7672 (idir==DIRECTION_K))
7674 force_zero_on_walls(rappel_ft[DIRECTION_K]);
7677 if (kmin + nz == nktot)
7690 for (
int k=nzdeb; k < nz; k++)
7691 for (
int j=0; j< ny; j++)
7692 for (
int i=0; i < nx; i++)
7696 const int icompo =
num_compo_[num_elem_zvdf];
7697 bool continuous =
false;
7698 for (
int c=0; c<nb_parts_continuous; c++)
7700 if (icompo == list_continuous_phase[c])
7706 if ((icompo==-1) || (continuous) )
7708 rappel_ft[idir](i,j,k) = 0.;
7712 int id_bulle = decodeur_num_compo[icompo];
7719 const int ibulle_reelle = decoder_numero_bulle(-ighost);
7722 id_bulle = ibulle_reelle;
7727 assert(id_bulle>=0);
7729 const double f = individual_forces(id_bulle,idir);
7730 integrated_forces(id_bulle,idir) += vol*f;
7731 integration_cells_per_bubble(id_bulle,idir) +=1;
7738 const double f = individual_forces(id_bulle,idir);
7739 rappel_ft[idir](i,j,k) = f;
7746 for (
int idir=0; idir < 3; idir++)
7750 individual_forces(ib, idir) = integrated_forces(ib, idir) / (vol*integration_cells_per_bubble(ib,idir));
7758 const IJK_Field_double& indic,
7759 const DoubleTab& individual_forces,
7760 const ArrOfDouble& volume_reel,
7761 const DoubleTab& position,
7762 const double coef_rayon_force_rappel)
7765 noms_forces.dimensionner_force(3);
7766 for (
int idir=0; idir < 3; idir++)
7768 noms_forces[idir] =
"";
7771 double xir = position(ib,0);
7772 double yir = position(ib,1);
7773 double zir = position(ib,2);
7774 double r2 = coef_rayon_force_rappel*coef_rayon_force_rappel*pow((volume_reel[ib]*3)/(4.*M_PI), 2./3.);
7775 noms_forces[idir] +=
"+((";
7776 noms_forces[idir] +=
Nom(
"(X-(")+
Nom(xir,
"%g")+
Nom(
"))*(X-(")+
Nom(xir,
"%g")+
Nom(
"))");
7777 noms_forces[idir] +=
Nom(
"+(Y-(")+
Nom(yir,
"%g")+
Nom(
"))*(Y-(")+
Nom(yir,
"%g")+
Nom(
"))");
7778 noms_forces[idir] +=
Nom(
"+(Z-(")+
Nom(zir,
"%g")+
Nom(
"))*(Z-(")+
Nom(zir,
"%g")+
Nom(
"))");
7779 noms_forces[idir] +=
Nom(
"-")+
Nom(r2,
"%g");
7783 noms_forces[idir] +=
")_LT_0.)*(1.-ff)*((1.-ff)_GT_0.000001)*(" ;
7784 noms_forces[idir] +=
Nom(individual_forces(ib,idir),
"%g");
7785 noms_forces[idir] +=
")" ;
7789 for (
int idir=0; idir < 3; idir++)
7797 const int ibulle_reelle = decoder_numero_bulle(-ighost);
7802 double r2 = coef_rayon_force_rappel*coef_rayon_force_rappel*pow((volume_reel[ibg]*3)/(4.*M_PI), 2./3.);
7803 noms_forces[idir] +=
"+((";
7804 noms_forces[idir] +=
Nom(
"(X-(")+
Nom(xir,
"%g")+
Nom(
"))*(X-(")+
Nom(xir,
"%g")+
Nom(
"))");
7805 noms_forces[idir] +=
Nom(
"+(Y-(")+
Nom(yir,
"%g")+
Nom(
"))*(Y-(")+
Nom(yir,
"%g")+
Nom(
"))");
7806 noms_forces[idir] +=
Nom(
"+(Z-(")+
Nom(zir,
"%g")+
Nom(
"))*(Z-(")+
Nom(zir,
"%g")+
Nom(
"))");
7807 noms_forces[idir] +=
Nom(
"-")+
Nom(r2,
"%g");
7810 noms_forces[idir] +=
")_LT_0.)*(1.-ff)*((1.-ff)_GT_0.000001)*(" ;
7811 noms_forces[idir] +=
Nom(individual_forces(ibulle_reelle,idir),
"%g");
7812 noms_forces[idir] +=
")" ;
7814 set_field_data(rappel[idir], noms_forces[idir], indic, 0. );
7819 const IJK_Field_double& indic,
7820 const ArrOfDouble& volume_reel,
7821 const DoubleTab& position)
const
7823 Nom nom_indicatrices_np;
7824 nom_indicatrices_np =
"";
7831 double r=pow((volume_reel[ib]*3)/(4.*M_PI), 1./3.);
7839 nom_indicatrices_np +=
"+((";
7840 nom_indicatrices_np +=
Nom(
"((X-(")+
Nom(x0,
"%g")+
Nom(
"))*(X-(")+
Nom(x0,
"%g")+
Nom(
"))/(")+
Nom(a,
"%g")+
Nom(
"*")+
Nom(a,
"%g")+
Nom(
"))");
7841 nom_indicatrices_np +=
Nom(
"+((Y-(")+
Nom(y0,
"%g")+
Nom(
"))*(Y-(")+
Nom(y0,
"%g")+
Nom(
"))/(")+
Nom(b,
"%g")+
Nom(
"*")+
Nom(b,
"%g")+
Nom(
"))");
7842 nom_indicatrices_np +=
Nom(
"+((Z-(")+
Nom(z0,
"%g")+
Nom(
"))*(Z-(")+
Nom(z0,
"%g")+
Nom(
"))/(")+
Nom(c,
"%g")+
Nom(
"*")+
Nom(c,
"%g")+
Nom(
"))");
7843 nom_indicatrices_np +=
Nom(
"-1.0");
7845 nom_indicatrices_np +=
")_LT_0.)*(1.0)" ;
7856 double r=pow((volume_reel[ibg]*3)/(4.*M_PI), 1./3.);
7864 nom_indicatrices_np +=
"+((";
7865 nom_indicatrices_np +=
Nom(
"((X-(")+
Nom(x0,
"%g")+
Nom(
"))*(X-(")+
Nom(x0,
"%g")+
Nom(
"))/(")+
Nom(a,
"%g")+
Nom(
"*")+
Nom(a,
"%g")+
Nom(
"))");
7866 nom_indicatrices_np +=
Nom(
"+((Y-(")+
Nom(y0,
"%g")+
Nom(
"))*(Y-(")+
Nom(y0,
"%g")+
Nom(
"))/(")+
Nom(b,
"%g")+
Nom(
"*")+
Nom(b,
"%g")+
Nom(
"))");
7867 nom_indicatrices_np +=
Nom(
"+((Z-(")+
Nom(z0,
"%g")+
Nom(
"))*(Z-(")+
Nom(z0,
"%g")+
Nom(
"))/(")+
Nom(c,
"%g")+
Nom(
"*")+
Nom(c,
"%g")+
Nom(
"))");
7868 nom_indicatrices_np +=
Nom(
"-1.0");
7870 nom_indicatrices_np +=
")_LT_0.)*(1.0)" ;
7873 Cerr <<
"Setting indicatrice_non_perturbe :" << nom_indicatrices_np << finl;
7875 set_field_data(indic_np, nom_indicatrices_np, indic, 0. );
7884 const DoubleTab& gravite,
7885 const double delta_rho,
7889 const bool parcourir
7898 statistics().create_custom_counter(
"Calcul Indicatrice Next",2,
"IJK");
7899 statistics().begin_count(
"Calcul Indicatrice Next",statistics().get_last_opened_counter_level()+1);
7910 statistics().end_count(
"Calcul Indicatrice Next");
7984 for (
int c=0; c < 3; c++)
8007 for (
int c=0; c<3; c++)
8019 for (
int bary_compo = 0; bary_compo < 3; bary_compo++)
8028 for (
int face_dir = 0; face_dir < 3; face_dir++)
8030 for (
int bary_compo = 0; bary_compo < 2; bary_compo++)
8040 for (
int d = 0; d < 3; d++)
8051 for (
int face_dir = 0; face_dir < 3; face_dir++)
8053 for (
int bary_compo = 0; bary_compo < 2; bary_compo++)
8059 statistics().end_count(
"Calcul Indicatrice Next");
8066 IJK_Field_double& indicatrice_intermediaire_ft,
8067 IJK_Field_double& indicatrice_intermediaire_ns,
8068 IJK_Field_vector3_double& indicatrice_surfacique_intermediaire_face_ft,
8069 IJK_Field_vector3_double& indicatrice_surfacique_intermediaire_face_ns,
8070 const bool parcourir
8076 statistics().create_custom_counter(
"Calcul Indicatrice Next",2,
"IJK");
8077 statistics().begin_count(
"Calcul Indicatrice Next",statistics().get_last_opened_counter_level()+1);
8081 indicatrice_intermediaire_ft.
data() = 1.;
8082 indicatrice_intermediaire_ns.
data() = 1.;
8088 statistics().end_count(
"Calcul Indicatrice Next");
8106 ref_ijk_ft_->eq_ns().redistrib_from_ft_elem().redistribute(
8107 indicatrice_intermediaire_ft, indicatrice_intermediaire_ns);
8114 ref_ijk_ft_->eq_ns().get_redistribute_from_splitting_ft_faces(
8115 indicatrice_surfacique_intermediaire_face_ft,
8116 indicatrice_surfacique_intermediaire_face_ns);
8119 statistics().end_count(
"Calcul Indicatrice Next");
8126 indicatrice_ft_test_.echange_espace_virtuel(indicatrice_ft_test_.ghost());
8131 groups_indicatrice_ft_test_.echange_espace_virtuel();
8133 SChaine group_indic;
8140 for (
int k = 0; k < nk; k++)
8141 for (
int j = 0; j < nj; j++)
8142 for (
int i = 0; i < ni; i++)
8147 indic <<
"(" << i <<
"," << j <<
"," << k <<
") : " <<
indicatrice_ft_[
next()](i, j, k) <<
" VS "
8148 << indicatrice_ft_test_(i, j, k) << finl;
8150 for (
int igroup = 0; igroup < nb_grps; igroup++)
8155 group_indic <<
"groupe = " << igroup <<
" (" << i <<
"," << j <<
"," << k
8157 << groups_indicatrice_ft_next_test_[igroup](i, j, k) << finl;
8165 Cerr <<
"Probleme_FTD_IJK_base:: calcul de l'indicatrice faux ! (iteration " << tstep_ <<
" sur " << nb_timesteps_ <<
" )"
8174 int current_time = next_time ?
next() :
old();
8175 int other_time = next_time ?
old() :
next();
8177 double current_indicatrice = next_time ?
In(i,j,k) :
I(i,j,k);
8178 double other_time_indicatrice = next_time ?
I(i,j,k) :
In(i,j,k);
8180 if ((
I(i,j,k) == 0.) && (
In(i,j,k) == 0.))
8184 else if ((
I(i,j,k) == 1.) && (
In(i,j,k) == 1.))
8188 else if (current_indicatrice == 0.)
8198 if (other_time_bary_compo > .5)
8208 else if (current_indicatrice == 1.)
8218 other_time_bary_compo =
opposing_barycentre(other_time_bary_compo, other_time_indicatrice);
8219 if (other_time_bary_compo > .5)
8237 assert(bary >= 0 && bary <= 1.);
8243 int current_time = next_time ?
next() :
old();
8244 int other_time = next_time ?
old() :
next();
8248 double current_indicatrice_surfacique = next_time ? next_indicatrice_surfacique : old_indicatrice_surfacique;
8249 double other_time_indicatrice_surfacique = next_time ? old_indicatrice_surfacique : next_indicatrice_surfacique;
8251 if (face_dir == bary_compo)
8260 int compo2D = (face_dir == 0) ? ((bary_compo == 2) ? 0 : ((bary_compo == 1) ? 1 : -1)) :
8261 ((face_dir == 1) ? ((bary_compo == 0) ? 0 : ((bary_compo == 2) ? 1 : -1)) :
8262 ((face_dir == 2) ? ((bary_compo == 1) ? 0 : ((bary_compo == 0) ? 1 : -1)) :
8264 assert(compo2D >= 0);
8265 if ((old_indicatrice_surfacique == 0.) && (next_indicatrice_surfacique == 0.))
8269 else if ((old_indicatrice_surfacique == 1.) && (next_indicatrice_surfacique == 1.))
8273 else if (current_indicatrice_surfacique == 0.)
8283 if (other_time_bary_compo > .5)
8293 else if (current_indicatrice_surfacique == 1.)
8303 other_time_bary_compo =
opposing_barycentre(other_time_bary_compo, other_time_indicatrice_surfacique);
8304 if (other_time_bary_compo > .5)
8322 assert(bary >= 0 && bary <= 1.);
8341 champs_compris_.switch_ft_fields();
8356 for (
int bary_compo = 0; bary_compo < 3; bary_compo++)
8365 for (
int face_dir = 0; face_dir < 3; face_dir++)
8367 for (
int bary_compo = 0; bary_compo < 2; bary_compo++)
8381 for (
int c = 0; c < 3; c++)
8397 const DoubleTab& gravite,
8398 const double delta_rho,
8406 ArrOfDouble potentiels_sommets;
8407 ArrOfDouble repulsions_sommets;
8410 potentiels_sommets, phi_par_compo);
8418 ArrOfDouble interfacial_source_term_sommet_dir, unite;
8425 for (
int dir = 0; dir < 3; dir++)
8428 interfacial_source_term_sommet_dir(som)=interfacial_source_term_sommet(som, dir);
8434 repulsions_sommets, repuls_par_compo);
8441 for (
int k = 0; k < nk; k++)
8442 for (
int j = 0; j < nj; j++)
8443 for (
int i = 0; i < ni; i++)
8451 ArrOfDouble& potentiels_sommets,
8452 ArrOfDouble& repulsions_sommets,
8453 const DoubleTab& gravite,
8454 const double delta_rho,
8469 potentiels_sommets = courbure;
8473 const double vol_cell = dxi * dxj * dxk;
8474 const DoubleTab& sommets = mesh.
sommets();
8475 const int nb_som = potentiels_sommets.
size_array();
8481 const ArrOfDouble& sigma_sommets =
maillage_ft_ijk_.Surfactant_facettes().get_sigma_sommets();
8482 for (
int i = 0; i < nb_som; i++)
8484 potentiels_sommets[i] *= sigma_sommets[i];
8489 potentiels_sommets *= sigma;
8500 ArrOfIntFT compo_connexe_sommets;
8505 DoubleTab deplacement;
8507 calculer_deplacement_from_code_compo_connexe_negatif(
8511 for (
int i = 0; i < nb_som; i++)
8513 double correction_potentiel_deplacement = 0.;
8514 if (compo_connexe_sommets[i] < 0)
8517 correction_potentiel_deplacement = -(
8518 gravite(0,0) * deplacement(i, 0) +
8519 gravite(0,1) * deplacement(i, 1) +
8520 gravite(0,2) * deplacement(i, 2)
8525 gravite(0,0) * coord[0] +
8526 gravite(0,1) * coord[1] +
8527 gravite(0,2) * coord[2] +
8528 correction_potentiel_deplacement
8530 potentiels_sommets[i] -= p * delta_rho;
8537 double vrx = DMAXFLOAT;
8538 double vry = DMAXFLOAT;
8539 double vrz = DMAXFLOAT;
8540 DoubleTab vr_to_closer;
8546 const int nznodes = z_grid_nodes.
size_array();
8547 const double zmin = z_grid_nodes[0];
8548 const double zmax = z_grid_nodes[nznodes - 1];
8549 double zsmin = (zmax - zmin) / 2.;
8550 for (
int i = 0; i < nb_som; i++)
8556 vrx = vr_to_closer(i, 0);
8557 vry = vr_to_closer(i, 1);
8558 vrz = vr_to_closer(i, 2);
8565 const double z = sommets(i, 2);
8566 double dzs = std::min(z - zmin, zmax - z);
8567 zsmin = std::min(zsmin, dzs);
8582 d = std::min(d, dzs);
8587 potentiels_sommets[i] -= phi;
8588 repulsions_sommets[i] = -phi * vol_cell;
8590 double dmin_local = dmin;
8596 if (fabs(dmin_local - dmin) < 1e-12)
8604 Cerr <<
"Warning. There were equalities ("
8606 <<
") in dmin computation. "
8607 <<
"2 equalities seems possible if it's between two markers "
8608 "separated by a proc frontieer. "
8623 envoyer_broadcast(vrx, iproc);
8624 envoyer_broadcast(vry, iproc);
8625 envoyer_broadcast(vrz, iproc);
8633 double norm_ur = std::sqrt(vrx * vrx + vry * vry + vrz * vrz);
8634 const double dt = ref_ijk_ft_->schema_temps_ijk().get_timestep();
8635 double cfl = dt * std::min(std::min(std::fabs(vrx) / dx, std::fabs(vry) / dy), std::fabs(vrz) / dz);
8636 int reset = (!
reprise_) && (itstep == 0);
8638 "_dmin.out",
"tstep\ttime\t\tdmin\tzsmin\tur\tCFL_local", reset
8640 fic << itstep <<
" " << time <<
" " << dmin <<
" " << zsmin;
8641 fic <<
" " << norm_ur <<
" " << cfl << finl;
8646 potentiels_sommets *= dxi * dxj * dxk;
8649 double lg_euler = std::max(std::max(dxi, dxj), dxk);
8651 const double sqrt_3 = 1.7320508075688772;
8652 bool ok = (lg_lagrange > lg_euler * sqrt_3);
8653 Cerr <<
"Test de la taille de maille eulerienne : lg_lagrange=" << lg_lagrange <<
" > (lg_euler=" << lg_euler
8654 <<
")*1.7320... ? " << (ok ?
"ok" :
"ko") <<
" (ratio : " << lg_lagrange / (lg_euler * sqrt_3)
8655 << (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)