TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Post_Processing_Hydrodynamic_Forces.cpp
1/****************************************************************************
2* Copyright (c) 2025, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
12* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
13*
14*****************************************************************************/
15
16#include <Post_Processing_Hydrodynamic_Forces.h>
17
18#include <Convection_Diffusion_Temperature_FT_Disc.h>
19#include <Transport_Interfaces_FT_Disc.h>
20#include <Navier_Stokes_FT_Disc.h>
21#include <Fluide_Incompressible.h>
22#include <Connex_components_FT.h>
23#include <Solid_Particle_base.h>
24#include <Fluide_Diphasique.h>
25#include <Maillage_FT_Disc.h>
26#include <Matrice_Dense.h>
27#include <Domaine_VDF.h>
28#include <TRUSTTab.h>
29#include <Objet_U_With_Params.h>
30
31
32Implemente_instanciable(Post_Processing_Hydrodynamic_Forces,
33 "Post_Processing_Hydrodynamic_Forces",Objet_U_With_Params);
34
35
36Sortie& Post_Processing_Hydrodynamic_Forces::printOn(Sortie& os) const
37{
38 Cerr << "Error: Post_Processing_Hydrodynamic_Forces::printOn is not implemented." << finl;
39 Process::exit();
40 return os;
41}
42
43void Post_Processing_Hydrodynamic_Forces::set_param(Param& p) const
44{
45 p.ajouter_flag("compute_hydrodynamic_forces", &is_compute_forces_);
46 p.ajouter_flag("compute_heat_transfer", &is_compute_heat_transfer_);
47 p.ajouter_flag("compute_stokes_theoretical_forces", &is_compute_stokes_theoretical_forces_);
48 p.ajouter_flag("post_process_pressure_fa7", &is_post_process_pressure_fa7_);
49 p.ajouter_flag("post_process_pressure_force_fa7", &is_post_process_pressure_force_fa7_);
50 p.ajouter_flag("post_process_friction_fa7", &is_post_process_friction_force_fa7_);
51 p.ajouter_flag("post_process_stress_tensor_fa7", &is_post_process_stress_tensor_fa7_);
52 p.ajouter("interpolation_distance_pressure_P1", &interpolation_distance_pressure_P1_);
53 p.ajouter("interpolation_distance_pressure_P2", &interpolation_distance_pressure_P2_);
54 p.ajouter("interpolation_distance_temperature_P1", &interpolation_distance_temperature_P1_);
55 p.ajouter("interpolation_distance_temperature_P2", &interpolation_distance_temperature_P2_);
56 p.ajouter("interpolation_distance_gradU_P1", &interpolation_distance_gradU_P1_);
57 p.ajouter("interpolation_distance_gradU_P2", &interpolation_distance_gradU_P2_);
58 p.ajouter_non_std("method_pressure_force_computation", (this), Param::REQUIRED);
59 p.ajouter_non_std("method_friction_force_computation", (this), Param::REQUIRED);
60 p.ajouter_non_std("location_stress_tensor", (this));
61}
62
63int Post_Processing_Hydrodynamic_Forces::lire_motcle_non_standard(const Motcle& mot, Entree& is)
64{
65
66 if (mot=="method_pressure_force_computation")
67 {
68 Motcles mots;
69 mots.add("trilinear_linear");
70 Motcle motbis;
71 is >> motbis;
72 Cerr << "Reading methode_calcul_force_pressure : " << motbis << finl;
73 const int r = mots.search(motbis);
74 switch(r)
75 {
76 case 0:
77 method_pressure_force_computation_ = Method_pressure_force_computation::TRILINEAR_LINEAR;
78 break;
79 default:
80 Cerr << "Error " << mots << "was expected whereas " << motbis <<" has been found."<< finl;
81 Process::exit();
82 }
83 return 1;
84 }
85 else if (mot=="method_friction_force_computation")
86 {
87 Motcles mots;
88 mots.add("trilinear_linear_complete_tensor");
89 mots.add("trilinear_linear_projected_tensor");
90 Motcle motbis;
91 is >> motbis;
92 Cerr << "Reading methode_calcul_force_frottements : " << motbis << finl;
93 const int r = mots.search(motbis);
94 switch(r)
95 {
96 case 0:
97 method_friction_force_computation_ =
98 Method_friction_force_computation::TRILINEAR_LINEAR_COMPLET_TENSOR;
99 break;
100 case 1:
101 method_friction_force_computation_ =
102 Method_friction_force_computation::TRILINEAR_LINEAR_PROJECTED_TENSOR;
103 break;
104 default:
105 Cerr << "Error " << mots << "was expected whereas " << motbis
106 << " has been found." << finl;
107 Process::exit();
108 }
109 return 1;
110 }
111 else if (mot=="location_stress_tensor")
112 {
113 Motcles mots;
114 mots.add("faces_normale_x");
115 mots.add("elements");
116 Motcle motbis;
117 is >> motbis;
118 Cerr << "Reading methode_interpolation : " << motbis << finl;
119 const int r = mots.search(motbis);
120 switch(r)
121 {
122 case 0:
123 location_stress_tensor_ = Location_stress_tensor::FACES_NORMALE_X;
124 break;
125 case 1:
126 location_stress_tensor_ = Location_stress_tensor::ELEMENTS;
127 break;
128 default:
129 Cerr << "Error " << mots << "was expected whereas " << motbis <<" has been found."<< finl;
130 Process::exit();
131 }
132 return 1;
133 }
134 else
135 {
136 Cerr << mot << " is not a keyword understood by " << que_suis_je() <<
137 " in lire_motcle_non_standard"<< finl;
138 Process::exit();
139 }
140 return -1;
141}
142
143
144void Post_Processing_Hydrodynamic_Forces::compute_hydrodynamic_forces()
145{
146 if (!flag_force_computation_)
147 {
148 Cerr << "Post_Processing_Hydrodynamic_Forces::compute_hydrodynamic_forces" << finl;
149 if (dimension!=3)
150 Process::exit("Post_Processing_Hydrodynamic_Forces::compute_hydrodynamic_forces"
151 " only implemented for 3D configurations.");
152 Transport_Interfaces_FT_Disc& eq_transport = ptr_eq_transport_.valeur();
153 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
154 const Maillage_FT_Disc& mesh = eq_transport.maillage_interface_pour_post();
155 const int nb_fa7 = mesh.nb_facettes();
156 IntVect compo_connexes_fa7(nb_fa7); // Init a zero
157 int n = search_connex_components_local_FT(mesh, compo_connexes_fa7);
158 int nb_particles_tot=compute_global_connex_components_FT(mesh, compo_connexes_fa7, n);
159
160 resize_and_init_tables(nb_particles_tot);
161
162 if (nb_fa7>0)
163 {
164 resize_data_fa7(nb_fa7);
165 resize_coord_neighbor_fluid_fa7(nb_fa7);
166
167 const ArrOfDouble& fa7_surface = mesh.get_update_surface_facettes();
168 const DoubleTab& tab_fa7_normal = mesh.get_update_normale_facettes();
169 const DoubleTab& gravity_center_fa7=mesh.get_gravity_center_fa7();
170
171 const OBS_PTR(Convection_Diffusion_Temperature_FT_Disc)& eq_thermal =
172 eq_ns.variables_internes().ref_equation_mpoint_;
173 bool is_discr_elem_diph=false;
174 Convection_Diffusion_Temperature_FT_Disc::Thermal_correction_discretization_method
175 thermal_correction_discretization_method=Convection_Diffusion_Temperature_FT_Disc::
176 Thermal_correction_discretization_method::P1;
177 if (eq_thermal)
178 {
179 thermal_correction_discretization_method=
180 eq_thermal.valeur().get_thermal_correction_discretization_method();
181 is_discr_elem_diph=thermal_correction_discretization_method==
182 Convection_Diffusion_Temperature_FT_Disc::
183 Thermal_correction_discretization_method::ELEM_DIPH;
184 }
185
186 if (is_discr_elem_diph)
187 list_elem_diph_.resize(nb_fa7,2); // EB (j,0) : num_elem_diph, (j,1) : num compo associee
188
189 const IntTab& particles_eulerian_id_number=eq_ns.get_particles_eulerian_id_number();
190 compute_neighbors_coordinates_fluid_fa7(nb_fa7,
191 is_discr_elem_diph,
192 gravity_center_fa7,
193 mesh,
194 tab_fa7_normal,
195 particles_eulerian_id_number);
196
197 // ----------------------------- Pressure force -----------------------------
198 compute_pressure_force_trilinear_linear(nb_fa7,mesh,thermal_correction_discretization_method,
199 compo_connexes_fa7, fa7_surface, tab_fa7_normal);
200 // ----------------------------- Friction force -----------------------------
201 if (method_friction_force_computation_==Method_friction_force_computation::
202 TRILINEAR_LINEAR_COMPLET_TENSOR)
203 {
204 compute_friction_force_complet_tensor(nb_fa7, mesh, compo_connexes_fa7,
205 fa7_surface, tab_fa7_normal);
206 }
207 else if (method_friction_force_computation_==Method_friction_force_computation::
208 TRILINEAR_LINEAR_PROJECTED_TENSOR) // see Butaye et al., Computers and Fluids, 2023.
209 {
210 compute_friction_force_projected_tensor(nb_fa7, mesh, compo_connexes_fa7,
211 fa7_surface, tab_fa7_normal);
212 }
213 }
214
215 mp_sum_for_each_item(total_pressure_force_);
216 mp_sum_for_each_item(total_friction_force_);
217 mp_sum_for_each_item(proportion_fa7_ok_UP2_);
218 mp_sum_for_each_item(total_surface_interf_);
219 mp_sum_for_each_item(prop_P2_fluid_compo_);
220 mp_sum_for_each_item(U_P2_moy_);
221 mp_sum_for_each_item(Nb_fa7_tot_par_compo_);
222
223 compute_U_P2_moy(nb_particles_tot);
224 compute_proportion_fa7_ok_and_is_fluid_P2(nb_particles_tot);
225
226 raise_the_flag();
227 }
228}
229
230int Post_Processing_Hydrodynamic_Forces::elem_faces_for_interp(int num_elem, int i) const
231{
232 const Domaine_VF& domain_vf=ref_cast(Domaine_VF,ptr_eq_ns_.valeur().domaine_dis());
233 if (num_elem<0 || num_elem>domain_vf.elem_faces().dimension_tot(0))
234 {
235 Journal() << "WARNING : IMPOSSIBLE to access the face of the element " << num_elem <<
236 " in the direction " << i << finl;
237 return -1;
238 }
239 return domain_vf.elem_faces(num_elem, i);
240}
241
242int Post_Processing_Hydrodynamic_Forces::face_voisins_for_interp(int num_face,int i) const
243{
244 const Domaine_VF& domain_vf=ref_cast(Domaine_VF,ptr_eq_ns_.valeur().domaine_dis());
245
246 if (num_face<0 || num_face > domain_vf.face_voisins().dimension_tot(0))
247 {
248 Journal() << "WARNING : IMPOSSIBLE to access the neighbor element of the face "
249 << num_face << " in the direction " << i << finl;
250 return -1;
251 }
252 return domain_vf.face_voisins(num_face,i);
253}
254
255int Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_elem(
256 const DoubleTab& valeurs_champ,
257 DoubleTab& coord,
258 DoubleTab& resu)
259{
260 return trilinear_interpolation_elem(valeurs_champ, coord, resu, 0,
261 Convection_Diffusion_Temperature_FT_Disc::Thermal_correction_discretization_method::P1);
262}
263
264// EB
265/*! @brief Interpolation trilineaire d'un champs "valeurs_champs" aux coordonnees coord a partir des 8 elements (4 en 2D) voisins les plus proches. Valeurs_chaamps contient
266 * typiquement le champ de presssion ou le champ de temperature.
267 * On utilise egalement cette fonction pour sauvegarder la liste des elments P1 ainsi que l'indiatrice des elements P2
268 */
269int Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_elem(
270 const DoubleTab& valeurs_champ,
271 DoubleTab& coord,
272 DoubleTab& resu,
273 const int is_P2,
274 const Convection_Diffusion_Temperature_FT_Disc::Thermal_correction_discretization_method
275 thermal_correction_discr_method)
276{
277
278 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, ptr_eq_ns_.valeur().domaine_dis());
279 int nb_voisins=8; // neighbouring elements
280 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
281 Transport_Interfaces_FT_Disc& eq_transport = ptr_eq_transport_.valeur();
282 const Maillage_FT_Disc& mesh = eq_transport.maillage_interface();
283 const DoubleTab& phase_indicator_function = eq_transport.get_indicatrice().valeurs();
284 const int nb_fa7=coord.dimension(0);
285
286 const IntTab& particles_eulerian_id_number=eq_ns.get_particles_eulerian_id_number();
287 const int sauv_list_P1 = (thermal_correction_discr_method==
288 Convection_Diffusion_Temperature_FT_Disc::
289 Thermal_correction_discretization_method::P1) ? 1 : 0;
290 if (sauv_list_P1)
291 {
292 list_elem_P1_.resize(nb_fa7,2); // In (j,0) we save element id, in (j,1) wa save the particle id
293 list_elem_P1_=-1;
294 }
295
296 // if P1_ALL, we discretize over all elements used for the interpolation in P1
297 const int sauv_list_P1_all = (thermal_correction_discr_method==
298 Convection_Diffusion_Temperature_FT_Disc::
299 Thermal_correction_discretization_method::P1_ALL) ? 1 : 0;
300 int nb_elem_P1_all=0;
301
302 if (sauv_list_P1_all)
303 list_elem_P1_all_.resize(0,2); // In (j,0) we save element id, in (j,1) wa save the particle id
304
305 const DoubleTab& gravity_center_fa7=mesh.get_gravity_center_fa7();
306
307 for (int fa7=0; fa7<nb_fa7; fa7++)
308 {
309 if (!mesh.facette_virtuelle(fa7))
310 {
311 DoubleVect coord_elem_interp(dimension);
312 for (int dim=0; dim<dimension; dim++)
313 coord_elem_interp(dim)=coord(fa7,dim);
314 IntVect neighboring_elements(nb_voisins);
315 if (is_P2)
316 {
317 phase_indicator_function_P2_(fa7)=find_neighboring_elements(
318 coord_elem_interp,
319 neighboring_elements);
320 }
321 else
322 {
323 find_neighboring_elements(coord_elem_interp,neighboring_elements,
324 sauv_list_P1,
325 fa7);
326 if (sauv_list_P1)
327 {
328 int elem_diph=domain_vdf.domaine().chercher_elements(
329 gravity_center_fa7(fa7,0),
330 gravity_center_fa7(fa7,1),
331 gravity_center_fa7(fa7,2));
332 int particle_eulerian_id_number = particles_eulerian_id_number(elem_diph);
333 list_elem_P1_(fa7,1)= particle_eulerian_id_number;
334 }
335 }
336 int access_elem=1;
337 for (int i=0; i<nb_voisins; i++)
338 {
339 if (neighboring_elements(i)<0)
340 access_elem= 0;
341 if (sauv_list_P1_all)
342 {
343 int elem_voisin=neighboring_elements(i);
344 int elem_diph=domain_vdf.domaine().chercher_elements(
345 gravity_center_fa7(fa7,0),
346 gravity_center_fa7(fa7,1),
347 gravity_center_fa7(fa7,2));
348 int particle_eulerian_id_number = particles_eulerian_id_number(elem_diph);
349 if (elem_voisin>=0 && phase_indicator_function(elem_voisin)==1)
350 {
351 list_elem_P1_all_.append_line(elem_voisin,particle_eulerian_id_number);
352 nb_elem_P1_all++;
353 }
354 }
355 }
356 if (access_elem)
357 {
358 DoubleVect delta_i(dimension);
359 delta_i(0) = fabs(domain_vdf.dist_elem(neighboring_elements(0),
360 neighboring_elements(1), 0));
361 delta_i(1) = fabs(domain_vdf.dist_elem(neighboring_elements(0),
362 neighboring_elements(2), 1));
363 delta_i(2) = fabs(domain_vdf.dist_elem(neighboring_elements(0),
364 neighboring_elements(4), 2));
365 DoubleVect coord_elem_0(dimension);
366 for (int dim=0; dim<dimension; dim++) coord_elem_0(dim)=domain_vdf.xp(
367 neighboring_elements(0),dim);
368
369 double xfact=fabs((coord_elem_interp(0)-coord_elem_0(0))/delta_i(0));
370 double yfact=fabs((coord_elem_interp(1)-coord_elem_0(1))/delta_i(1));
371 double zfact=fabs((coord_elem_interp(2)-coord_elem_0(2))/delta_i(2));
372
373 resu(fa7)=(1-zfact)*(
374 (1-yfact)*(
375 (1-xfact)*(valeurs_champ(neighboring_elements(0))) +
376 xfact*(valeurs_champ(neighboring_elements(1)))
377 ) +
378 yfact*(
379 (1-xfact)*(valeurs_champ(neighboring_elements(2))) +
380 xfact*(valeurs_champ(neighboring_elements(3)))
381 )
382 ) +
383 zfact*(
384 (1-yfact)*(
385 (1-xfact)*(valeurs_champ(neighboring_elements(4))) +
386 xfact*(valeurs_champ(neighboring_elements(5)))
387 ) +
388 yfact*(
389 (1-xfact)*(valeurs_champ(neighboring_elements(6))) +
390 xfact*(valeurs_champ(neighboring_elements(7)))
391 )
392 );
393 }
394 else
395 resu(fa7)=-1e15;
396 }
397 }
398 if (sauv_list_P1_all)
399 {
400 list_elem_P1_all_.resize(nb_elem_P1_all,2);
401 }
402 return 1;
403}
404
405
406int Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_face(
407 const DoubleTab& valeurs_champ,
408 DoubleTab& coord,
409 DoubleTab& resu)
410{
411 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
412 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, eq_ns.domaine_dis());
413 const Transport_Interfaces_FT_Disc& eq_transport = ptr_eq_transport_.valeur();
414 const Maillage_FT_Disc& mesh = eq_transport.maillage_interface();
415
416 int nb_neighbors=8;
417 for (int fa7=0; fa7<coord.dimension(0); fa7++)
418 {
419 if (!mesh.facette_virtuelle(fa7))
420 {
421 // On recupere les faces voisines : les faces euleriennes qui vont servir a l'interpolation
422
423 DoubleVect coord_elem_interp(dimension); // coordinate of the fa7 gravity center
424 for (int dim=0; dim<dimension; dim++)
425 coord_elem_interp(dim)=coord(fa7,dim);
426 IntTab faces_voisines(dimension,nb_neighbors); // 8 elements voisins au point de coordonnees coord. L'element elem est inclu dedans.
427 find_neighboring_faces_xyz(coord_elem_interp,faces_voisines);
428 int acces_faces=1;
429
430 for (int dim=0; dim<dimension; dim++)
431 {
432 for (int i=0; i<nb_neighbors; i++)
433 {
434 if (faces_voisines(dim,i)<0)
435 acces_faces= 0; // s'il y a eu un bug dans la recuperation des faces (pbm d'acces : zone de joint, bord), on renvoie 0 et la fa7 n'est pas prise en compte dans le calcul
436 }
437 }
438
439 if (acces_faces)
440 {
441 DoubleVect xfact(dimension);
442 DoubleVect yfact(dimension);
443 DoubleVect zfact(dimension);
444 DoubleTab Delta_(dimension);
445 Delta_(0)=fabs(domain_vdf.dist_face(faces_voisines(0,0),faces_voisines(0,1),0));
446 Delta_(1)=fabs(domain_vdf.dist_face(faces_voisines(0,0),faces_voisines(0,2),1));
447 Delta_(2)=fabs(domain_vdf.dist_face(faces_voisines(0,0),faces_voisines(0,4),2));
448
449 for (int dim=0; dim<dimension; dim++)
450 {
451 xfact(dim)=fabs((coord(fa7,0)-domain_vdf.xv(faces_voisines(dim,0),0))/Delta_(0));
452 yfact(dim)=fabs((coord(fa7,1)-domain_vdf.xv(faces_voisines(dim,0),1))/Delta_(1));
453 zfact(dim)=fabs((coord(fa7,2)-domain_vdf.xv(faces_voisines(dim,0),2))/Delta_(2));
454 }
455
456 for (int dim=0; dim<dimension; dim++)
457 {
458 if (xfact(dim)>1)
459 Cerr << "xfact > 1 pour le point " <<coord(fa7,0)<< " " << coord(fa7,1)<< " "
460 << coord(fa7,2) <<finl;
461 if (yfact(dim)>1)
462 Cerr << "yfact > 1 pour le point " <<coord(fa7,0)<< " " << coord(fa7,1)<< " "
463 << coord(fa7,2) <<finl;
464 if (zfact(dim)>1)
465 Cerr << "zfact > 1 pour le point " <<coord(fa7,0)<< " " << coord(fa7,1)<< " "
466 << coord(fa7,2) <<finl;
467 }
468
469 for (int dim=0; dim<dimension; dim++)
470 {
471 resu(fa7,dim)=(1.-zfact(dim))*(
472 (1.-yfact(dim))*(
473 (1.-xfact(dim))*(valeurs_champ(faces_voisines(dim,0))) +
474 xfact(dim)*(valeurs_champ(faces_voisines(dim,1)))
475 ) +
476 yfact(dim)*(
477 (1.-xfact(dim))*(valeurs_champ(faces_voisines(dim,2))) +
478 xfact(dim)*(valeurs_champ(faces_voisines(dim,3))))
479 ) +
480 zfact(dim)*(
481 (1.-yfact(dim))*(
482 (1.-xfact(dim))*(valeurs_champ(faces_voisines(dim,4))) +
483 xfact(dim)*(valeurs_champ(faces_voisines(dim,5)))
484 ) +
485 yfact(dim)*(
486 (1.-xfact(dim))*(valeurs_champ(faces_voisines(dim,6))) +
487 xfact(dim)*(valeurs_champ(faces_voisines(dim,7))))
488 );
489 }
490 }
491 else
492 {
493 for (int dim=0; dim<dimension; dim++)
494 resu(fa7,dim)=-1e15; // de cette maniere, on ne calcule pas la force pour la fa7 pour laquelle on n'a pas acces a P2
495 }
496 }
497 }
498 return 1;
499}
500
501int Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_gradU_face(
502 const DoubleTab& valeurs_champ,
503 DoubleTab& coord,
504 DoubleTab& resu)
505{
506 const Transport_Interfaces_FT_Disc& eq_transport = ptr_eq_transport_.valeur();
507 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
508 const Maillage_FT_Disc& mesh = eq_transport.maillage_interface();
509 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, eq_ns.domaine_dis());
510 const Domaine& domain = domain_vdf.domaine();
511 const DoubleTab& xv = domain_vdf.xv();
512 IntVect faces_elem_interp(2*dimension);
513 int nb_neighbors=8;
514
515 for (int fa7=0; fa7<coord.dimension(0); fa7++)
516 {
517 if (!mesh.facette_virtuelle(fa7))
518 {
519 const int elem=domain.chercher_elements(coord(fa7,0), coord(fa7,1), coord(fa7,2));
520 for (int dim=0; dim<dimension; dim++)
521 {
522 faces_elem_interp(dim)=elem_faces_for_interp(elem,dim);
523 faces_elem_interp(dimension+dim)=elem_faces_for_interp(elem,dimension+dim);
524 if (faces_elem_interp(dim)<0 || faces_elem_interp(dimension+dim)<0)
525 return 0;
526 }
527
528 DoubleTab coord_face(2*dimension,dimension);
529 for (int i=0; i<2*dimension; i++)
530 {
531 for (int dim=0; dim<dimension; dim++)
532 {
533 coord_face(i,dim)=domain_vdf.xv(faces_elem_interp(i),dim);
534 }
535 }
536
537 DoubleVect coord_elem_interp(dimension);
538 for (int dim=0; dim<dimension; dim++) coord_elem_interp(dim)=coord(fa7,dim);
539 IntVect faces_voisines(nb_neighbors); // 8 elements voisins au point de coordonnees coord. L'element elem est inclu dedans.
540 find_neighboring_faces(coord_elem_interp,faces_voisines,0);
541
542 for (int i=0; i<nb_neighbors; i++)
543 {
544 if (faces_voisines(i)<0)
545 return 0;
546 }
547
548 DoubleTab gradUx(nb_neighbors, dimension, dimension); // le x signifie que l'on calcule le tenseur en une facette dont la composante de la vitesse est en x
549 for (int face=0; face<8; face++) //on calcule le tenseur gradient de la vitesse pour chaque face voisine
550 {
551 // SCHEMA EN 2D
552 // ---- --- --- --- --- --- --- --- --- --- -
553 // | | | | | | | | | | |
554 // 3
555 // | | | | | | | | | | |
556 // ---- --- --- --- -7- -8- --- --- --- --- -
557 // | | | | | | | | | | |
558 // 1 x 2
559 // | | | | | | | | | | |
560 // ---- --- --- --- -5- -6- --- --- --- --- - y
561 // | | | | | | | | | | | ^
562 // 4 |
563 // | | | | | | | | | | | |
564 // ---- --- --- --- --- --- --- --- --- --- - o----->x
565 // z
566 // Evaluation du gradient sur les faces de normale x (choix arbitraire mais peut tuer une certaine symetrie)
567 // Les termes + designent les faces juxtaposees selon +z
568 // Les termes - designent les faces juxtaposees selon -z
569 // Les termes 57 et 86 pour w designent les faces de normale z secantes aux faces 5,7 et 8,6
570 // -- --
571 // | u_2-u_1 u_3-u_4 u_x+ - u_x- |
572 // | ---------- ------------ --------------- |
573 // | (x_2-x1) (y_3-y_4) (z_x+ - z_x-) |
574 // | |
575 // grad U = | 1 v_7-v_8 v_5 - v_6 1 v_5-v_7 v_6-v_8 1 v_7+ - v_7- v_8+ - v_8- v_5+ - v_5- v_6+ - v_6- |
576 // | - ( --------- + --------- ) - ( ------ + -------) - ( --------------- + --------------- + --------------- + --------------- ) |
577 // | 2 x_7-x_8 x_5 - x_6 2 y_5-y_7 y_6-y_8 4 (z_7+ - z_7-) (z_8+ - z_8-) (z_5+ - z_8-) (z_6+ - z_6-) |
578 // | |
579 // | 1 w57+ - w86+ w57- - w86- 1 w27-w29 w28-w30 w23-w25 w24-w26 1 w_57+ - w_57- w_86+ - w_86- |
580 // | - ( ----------- + ----------- ) - ( -------- + -------- + --------- + --------- ) - ( ------------ + --------------) |
581 // | 2 x57+ - x86+ x57- - x86- 4 y27-y29 y28-y30 y23-y25 y24-y26 2 z_57+ - z_57- z_86+ - z_86- |
582 // -- --
583 //
584 // Pour chaque face, il faut donc 30 faces voisines : 1-8, x+, x-, 5+, 5-, 6+, 6-, 7+, 7-, 8+, 8-, 57+, 57-, 86+, 86-
585 // Pour plus de facilites, on numerote : x+:9, x-:10, 5+:11, 5-:12, 6+:13, 6-:14, 7+:15, 7-:16, 8+:17, 8-:18, 57+:19, 57-:20, 86+:21, 86-:22
586 // 23-30 : les faces de normales z autour de x. Numerotation de bas en haut, de devant a derriere, de gauche a droite
587
588 int nb_voisinsx=30;
589 IntVect voisinsx(nb_voisinsx);
590 int num_facex=faces_voisines(face);
591 int elem_gauche=face_voisins_for_interp(num_facex, 0);
592 int elem_droite=face_voisins_for_interp(num_facex, 1);
593 int elem_haut_gauche=face_voisins_for_interp(elem_faces_for_interp(elem_gauche,2+dimension), 1);
594 int elem_haut_droite=face_voisins_for_interp(elem_faces_for_interp(elem_droite,2+dimension), 1);
595 int elem_bas_gauche=face_voisins_for_interp(elem_faces_for_interp(elem_gauche,2), 0);
596 int elem_bas_droite=face_voisins_for_interp(elem_faces_for_interp(elem_droite,2), 0);
597 int elem_avant_gauche=face_voisins_for_interp(elem_faces_for_interp(elem_gauche,1+dimension),1);
598 int elem_avant_droite=face_voisins_for_interp(elem_faces_for_interp(elem_droite,1+dimension),1);
599 int elem_arriere_gauche=face_voisins_for_interp(elem_faces_for_interp(elem_gauche,1),0);
600 int elem_arriere_droite=face_voisins_for_interp(elem_faces_for_interp(elem_droite,1),0);
601
602 voisinsx(0) = elem_faces_for_interp(elem_gauche,0);
603 voisinsx(1) = elem_faces_for_interp(elem_droite,0+dimension);
604 voisinsx(2) = elem_faces_for_interp(elem_avant_gauche,0+dimension);
605 voisinsx(3) = elem_faces_for_interp(elem_arriere_gauche,0+dimension);
606 voisinsx(4) = elem_faces_for_interp(elem_gauche,1);
607 voisinsx(5) = elem_faces_for_interp(elem_droite,1);
608 voisinsx(6) = elem_faces_for_interp(elem_gauche,1+dimension);
609 voisinsx(7) = elem_faces_for_interp(elem_droite,1+dimension);
610 voisinsx(8) = elem_faces_for_interp(elem_haut_gauche,0+dimension);
611 voisinsx(9) = elem_faces_for_interp(elem_bas_gauche,0+dimension);
612 voisinsx(10) = elem_faces_for_interp(elem_haut_gauche,1);
613 voisinsx(11) = elem_faces_for_interp(elem_bas_gauche,1);
614 voisinsx(12) = elem_faces_for_interp(elem_haut_droite,1);
615 voisinsx(13) = elem_faces_for_interp(elem_bas_droite,1);
616 voisinsx(14) = elem_faces_for_interp(elem_haut_gauche,1+dimension);
617 voisinsx(15) = elem_faces_for_interp(elem_bas_gauche,1+dimension);
618 voisinsx(16) = elem_faces_for_interp(elem_haut_droite,1+dimension);
619 voisinsx(17) = elem_faces_for_interp(elem_bas_droite,1+dimension);
620 voisinsx(18) = elem_faces_for_interp(elem_gauche,2+dimension);
621 voisinsx(19) = elem_faces_for_interp(elem_gauche,2);
622 voisinsx(20) = elem_faces_for_interp(elem_droite,2+dimension);
623 voisinsx(21) = elem_faces_for_interp(elem_droite,2);
624 voisinsx(22) = elem_faces_for_interp(elem_arriere_gauche,2);
625 voisinsx(23) = elem_faces_for_interp(elem_arriere_droite,2);
626 voisinsx(24) = elem_faces_for_interp(elem_avant_gauche,2);
627 voisinsx(25) = elem_faces_for_interp(elem_avant_droite,2);
628 voisinsx(26) = elem_faces_for_interp(elem_arriere_gauche,2+dimension);
629 voisinsx(27) = elem_faces_for_interp(elem_arriere_droite,2+dimension);
630 voisinsx(28) = elem_faces_for_interp(elem_avant_gauche,2+dimension);
631 voisinsx(29) = elem_faces_for_interp(elem_avant_droite,2+dimension);
632
633 for (int i=0; i<nb_voisinsx; i++)
634 {
635 if (voisinsx(i)<0)
636 return 0;
637 }
638
639 gradUx(face,0,0) = (valeurs_champ(voisinsx(1)) - valeurs_champ(voisinsx(0))) /
640 (xv(voisinsx(1),0) - xv(voisinsx(0),0));
641
642 gradUx(face,0,1) = (valeurs_champ(voisinsx(2)) - valeurs_champ(voisinsx(3))) /
643 (xv(voisinsx(2),1) - xv(voisinsx(3),1));
644
645 gradUx(face,0,2) = (valeurs_champ(voisinsx(8)) - valeurs_champ(voisinsx(9))) /
646 (xv(voisinsx(8),2) - xv(voisinsx(9),2));
647
648 gradUx(face,1,0) = 1./2.*((valeurs_champ(voisinsx(6)) - valeurs_champ(voisinsx(7))) /
649 (xv(voisinsx(6),0) - xv(voisinsx(7),0)) + (valeurs_champ(voisinsx(4)) -
650 valeurs_champ(voisinsx(5))) / (xv(voisinsx(4),0) - xv(voisinsx(5),0)));
651
652 gradUx(face,1,1) = 1./2.*((valeurs_champ(voisinsx(4)) - valeurs_champ(voisinsx(6))) /
653 (xv(voisinsx(4),1) - xv(voisinsx(6),1)) + (valeurs_champ(voisinsx(5)) -
654 valeurs_champ(voisinsx(7))) / (xv(voisinsx(5),1) - xv(voisinsx(7),1)));
655
656 gradUx(face,1,2) = 1./4.*((valeurs_champ(voisinsx(14)) - valeurs_champ(voisinsx(15))) /
657 (xv(voisinsx(14),2) - xv(voisinsx(15),2)) + (valeurs_champ(voisinsx(16)) -
658 valeurs_champ(voisinsx(17))) / (xv(voisinsx(16),2) - xv(voisinsx(17),2)) +
659 (valeurs_champ(voisinsx(10)) - valeurs_champ(voisinsx(11))) / (xv(voisinsx(10),2)
660 - xv(voisinsx(11),2)) + (valeurs_champ(voisinsx(12)) -
661 valeurs_champ(voisinsx(13))) / (xv(voisinsx(12),2) - xv(voisinsx(13),2)));
662
663 gradUx(face,2,0) = 1./2.*((valeurs_champ(voisinsx(18)) - valeurs_champ(voisinsx(20))) /
664 (xv(voisinsx(18),0) - xv(voisinsx(20),0)) + (valeurs_champ(voisinsx(19)) -
665 valeurs_champ(voisinsx(21))) / (xv(voisinsx(19),0) - xv(voisinsx(21),0)));
666
667 gradUx(face,2,1) = 1./4.*((valeurs_champ(voisinsx(26)) - valeurs_champ(voisinsx(28))) /
668 (xv(voisinsx(26),1) - xv(voisinsx(28),1)) + (valeurs_champ(voisinsx(27)) -
669 valeurs_champ(voisinsx(29))) / (xv(voisinsx(27),1) - xv(voisinsx(29),1)) +
670 (valeurs_champ(voisinsx(22)) - valeurs_champ(voisinsx(24))) / (xv(voisinsx(22),1) -
671 xv(voisinsx(24),1)) + (valeurs_champ(voisinsx(23)) - valeurs_champ(voisinsx(25))) /
672 (xv(voisinsx(23),1) - xv(voisinsx(25),1)));
673
674 gradUx(face,2,2) = 1./2.*((valeurs_champ(voisinsx(18)) - valeurs_champ(voisinsx(19))) /
675 (xv(voisinsx(18),2) - xv(voisinsx(19),2)) + (valeurs_champ(voisinsx(20)) -
676 valeurs_champ(voisinsx(21))) / (xv(voisinsx(20),2) - xv(voisinsx(21),2)));
677 }
678
679 double xfact;
680 double yfact;
681 double zfact;
682
683 DoubleTab Delta_x(dimension);
684
685 Delta_x(0)=fabs(domain_vdf.dist_face(faces_voisines(0),faces_voisines(1),0));
686 Delta_x(1)=fabs(domain_vdf.dist_face(faces_voisines(0),faces_voisines(2),1));
687 Delta_x(2)=fabs(domain_vdf.dist_face(faces_voisines(0),faces_voisines(4),2));
688
689 xfact=fabs((coord(fa7,0)-coord_face(0,0))/Delta_x(0));
690 yfact=fabs((coord(fa7,1)-coord_face(0,1))/Delta_x(1));
691 zfact=fabs((coord(fa7,2)-coord_face(0,2))/Delta_x(2));
692
693 for (int i=0; i<dimension; i++)
694 {
695 for (int j=0; j<dimension; j++)
696 {
697 resu(fa7,i,j)=(1-zfact)*(
698 (1-yfact)*( (1-xfact)*(gradUx(0,i,j)) + xfact*(gradUx(1,i,j)) ) +
699 yfact*( (1-xfact)*(gradUx(2,i,j)) + xfact*(gradUx(3,i,j)) )
700 ) +
701 zfact*(
702 (1-yfact)*( (1-xfact)*(gradUx(4,i,j)) + xfact*(gradUx(5,i,j)) ) +
703 yfact*( (1-xfact)*(gradUx(6,i,j)) + xfact*(gradUx(7,i,j))) );
704 }
705 }
706 }
707 }
708 return 1;
709}
710
711int Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_gradU_elem_P1(
712 const DoubleTab& valeurs_champ,
713 DoubleTab& coord,
714 DoubleTab& resu)
715{
716 Transport_Interfaces_FT_Disc& eq_transport = ptr_eq_transport_.valeur();
717 const DoubleTab& phase_indicator_function = eq_transport.get_indicatrice().valeurs();
718 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
719 const IntTab& particles_eulerian_id_number=eq_ns.get_particles_eulerian_id_number();
720 const Fluide_Diphasique& two_phase_fluid = eq_ns.fluide_diphasique();
721 const Maillage_FT_Disc& mesh = eq_transport.maillage_interface();
722 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, eq_ns.domaine_dis());
723 const DoubleTab& xv = domain_vdf.xv();
724 int nb_voisins=8;
725 const int id_fluid_phase=two_phase_fluid.get_id_fluid_phase();
726 const double mu_f =two_phase_fluid.fluide_phase(id_fluid_phase).
727 viscosite_dynamique().valeurs()(0, 0);
728 const double mu_p = two_phase_fluid.fluide_phase(1-id_fluid_phase).
729 viscosite_dynamique().valeurs()(0, 0);
730
731 for (int fa7=0; fa7<coord.dimension(0); fa7++)
732 {
733 if (!mesh.facette_virtuelle(fa7))
734 {
735 // find the 8 neighboring elements where the stress tensor will be computed
736 DoubleVect coord_elem_interp(dimension);
737 for (int dim=0; dim<dimension; dim++) coord_elem_interp(dim)=coord(fa7,dim);
738 IntVect elem_voisins(nb_voisins);
739 find_neighboring_elements(coord_elem_interp,elem_voisins);
740
741 for (int i=0; i<nb_voisins; i++)
742 {
743 if (elem_voisins(i)<0)
744 return 0;
745 }
746
747 // compute distance between neighboring meshes in all direction to compute
748 // interpolation coefficients
749 DoubleVect coord_elem_0(dimension);
750 DoubleVect delta_i(dimension);
751 delta_i(0) = fabs(domain_vdf.dist_elem(elem_voisins(0), elem_voisins(1), 0));
752 delta_i(1) = fabs(domain_vdf.dist_elem(elem_voisins(0), elem_voisins(3), 1));
753 delta_i(2) = fabs(domain_vdf.dist_elem(elem_voisins(0), elem_voisins(5), 2));
754
755 for (int dim=0; dim<dimension; dim++)
756 coord_elem_0(dim)=domain_vdf.xp(elem_voisins(0),dim);
757
758 double xfact=fabs((coord_elem_interp(0)-coord_elem_0(0))/delta_i(0));
759 double yfact=fabs((coord_elem_interp(1)-coord_elem_0(1))/delta_i(1));
760 double zfact=fabs((coord_elem_interp(2)-coord_elem_0(2))/delta_i(2));
761
762 // compute the velocity gradient tensor for each neighboring element
763 DoubleTab gradU(nb_voisins, dimension, dimension);
764 for (int elem=0; elem<nb_voisins; elem++)
765 {
766
767 // SCHEMA EN 2D
768 // FRONT ---- ---- --- --- --- --- --- --- --- FRONT
769 // | | | | | | | | | |
770 // 9 10
771 // | | | | | | | | | |
772 // ---- --- --- -5- -3- -6- --- --- ---
773 // | | | | | | | | | |
774 // LEFT 1 x 2 RIGHT
775 // | | | | | | | | | |
776 // ---- --- --- -7- -4- -8- --- --- --- y
777 // | | | | | | | | | | ^
778 // 11 12 |
779 // | | | | | | | | | | |
780 // BACK ---- --- --- --- --- --- --- --- --- --- BACK o----->x
781 // z
782 // Evaluation du gradient sur les faces de normale x (choix arbitraire mais peut tuer une certaine symetrie)
783 // Les termes + designent les faces juxtaposees selon +z
784 // Les termes - designent les faces juxtaposees selon -z
785 // Les termes 57 et 86 pour w designent les faces de normale z secantes aux faces 5,7 et 8,6
786 //
787 // ON CALCULE mu*gradU et pas gradU
788 // Pour plus de lisibilite, mu n'est pas reecrit pour chaque composante (mais sans ce mu_h local, la decomposition n'as plus d'interet)
789 // ON UTILISE LE MU HARMONIQUE LOCAL, IE : LE MU_HARMONIQUE AUX ARETES
790 //
791 // -- --
792 // | u_2-u_1 1 1 u_9-u_1 u_1-u_11 1 u_10-u_2 u_2-u_12 1 1 u_1+ - u_1 u_1 - u_1- 1 u_2+ - u_2 u_2 - u_2- |
793 // | ---------- -( - (--------- + --------)+ - (--------- + ---------) ) - ( - ( ----------- + ---------- ) + - ( ----------- + ---------- ) ) |
794 // | x_2-x1 2 2 y_9-y_1 y_1-y_11 2 y_10-y_2 y_2-y_12 2 2 z_1+ - z_1 z_1 - z_1- 2 z_2+ - z_2 z_2 - z_2- |
795 // | |
796 // grad U = | 1 1 v_5-v_3 v_3 - v_6 1 v_7-v_4 v_4-v_8 v_3-v_4 1 1 v_3+ - v_3 v_3 - v_3- 1 v_4+ - v_4 v_4 - v_4- |
797 // | - ( - ( -------- + ---------)+ - (-------- + --------) ) ------- - ( - ( ----------- + ---------- ) + - ( ----------- + ---------- ) ) |
798 // | 2 2 x_5-x_3 x_3 - x_6 2 x_7-x_4 x_4-x_8 y_3-y_4 2 2 z_3+ - z_3 z_3 - z_3- 2 z_4+ - z_4 z_4 - z_4- |
799 // | |
800 // | 1 1 w57+ - w34+ w34+ - w86+ 1 w57- - w34- w34- - w86- 1 1 w_1112+ - w_34+ w_34+ - w_910+ 1 w_1112- - w_34- w_34- - w_910- w_34+ - w_34- |
801 // | - ( - (---------- + ----------- ) + -(---------- + ------------) ) - ( - ( --------------- + -------------- ) + - ( --------------- + -------------- ) ) ------------- |
802 // | 2 2 x57+ - x34+ x34+ - x86+ 2 x57- - x34- x34- - x86- 2 2 y_1112+ - y_34+ y_34+ - y_910+ 2 y_1112- - y_34- y_34- - y_910- z_34+ - z_34- |
803 // -- --
804 //
805 // Pour chaque face, il faut donc 30 faces voisines : 1-12, 1+, 1-, 2+, 2-, 3+, 3-, 4+, 4-, 34+, 34-, 57+, 57-, 86+, 86-, 1112+, 1112-, 910+, 910-
806 // Pour plus de facilites, on numerote : 1+:13, 1-:14, 2+:15, 2-:16, 3+:17, 3-:18, 4+:19, 4-:20, 34+:21, 34-:22, 57+:23, 57-:24, 86+:25, 86-:26, 1112+:27, 1112-:28, 910+:29, 910-:30
807 // On identifie les faces voisines pour le calcul des composantes du tenseur
808
809 int nb_faces_voisines=30;
810 IntVect les_faces_voisines(nb_faces_voisines);
811 int elem_=elem_voisins(elem);
812 int elem_gauche = face_voisins_for_interp(elem_faces_for_interp(elem_,0),0);
813 int elem_droite = face_voisins_for_interp(elem_faces_for_interp(elem_,0+dimension),1);
814 int elem_haut=face_voisins_for_interp(elem_faces_for_interp(elem_,2+dimension), 1);
815 int elem_bas=face_voisins_for_interp(elem_faces_for_interp(elem_,2), 0);
816 int elem_avant=face_voisins_for_interp(elem_faces_for_interp(elem_,1+dimension),1);
817 int elem_arriere=face_voisins_for_interp(elem_faces_for_interp(elem_,1),0);
818
819 les_faces_voisines(0)=elem_faces_for_interp(elem_,0);
820 les_faces_voisines(1)=elem_faces_for_interp(elem_,0+dimension);
821 les_faces_voisines(2)=elem_faces_for_interp(elem_,1);
822 les_faces_voisines(3)=elem_faces_for_interp(elem_,1+dimension);
823 les_faces_voisines(4)=elem_faces_for_interp(elem_gauche,1+dimension);
824 les_faces_voisines(5)=elem_faces_for_interp(elem_droite,1+dimension);
825 les_faces_voisines(6)=elem_faces_for_interp(elem_gauche,1);
826 les_faces_voisines(7)=elem_faces_for_interp(elem_droite,1);
827 les_faces_voisines(8)=elem_faces_for_interp(elem_avant,0);
828 les_faces_voisines(9)=elem_faces_for_interp(elem_avant,0+dimension);
829 les_faces_voisines(10)=elem_faces_for_interp(elem_arriere,0);
830 les_faces_voisines(11)=elem_faces_for_interp(elem_arriere,0+dimension);
831 les_faces_voisines(12)=elem_faces_for_interp(elem_haut,0);
832 les_faces_voisines(13)=elem_faces_for_interp(elem_bas,0);
833 les_faces_voisines(14)=elem_faces_for_interp(elem_haut,0+dimension);
834 les_faces_voisines(15)=elem_faces_for_interp(elem_bas,0+dimension);
835 les_faces_voisines(16)=elem_faces_for_interp(elem_haut,1+dimension);
836 les_faces_voisines(17)=elem_faces_for_interp(elem_bas,1+dimension);
837 les_faces_voisines(18)=elem_faces_for_interp(elem_haut,1);
838 les_faces_voisines(19)=elem_faces_for_interp(elem_bas,1);
839 les_faces_voisines(20)=elem_faces_for_interp(elem_,2);
840 les_faces_voisines(21)=elem_faces_for_interp(elem_,2+dimension);
841 les_faces_voisines(22)=elem_faces_for_interp(elem_gauche,2);
842 les_faces_voisines(23)=elem_faces_for_interp(elem_gauche,2+dimension);
843 les_faces_voisines(24)=elem_faces_for_interp(elem_droite,2);
844 les_faces_voisines(25)=elem_faces_for_interp(elem_droite,2+dimension);
845 les_faces_voisines(26)=elem_faces_for_interp(elem_arriere,2);
846 les_faces_voisines(27)=elem_faces_for_interp(elem_arriere,2+dimension);
847 les_faces_voisines(28)=elem_faces_for_interp(elem_avant,2);
848 les_faces_voisines(29)=elem_faces_for_interp(elem_avant,2+dimension);
849
850 for (int i=0; i<nb_faces_voisines; i++)
851 {
852 if (les_faces_voisines(i)<0)
853 return 0;
854 }
855
856 int particle_id= particles_eulerian_id_number(elem);
857
858 gradU(elem,0,0) = ( (mu_p*mu_f/(mu_f-phase_indicator_function(elem_voisins(elem))*
859 (mu_f-mu_p)))*(valeurs_champ(les_faces_voisines(1)) -
860 valeurs_champ(les_faces_voisines(0))) /
861 (xv(les_faces_voisines(1),0) - xv(les_faces_voisines(0),0)));
862
863 gradU(elem,0,1) = 1./2.* ( 1./2.*( compute_viscosity_edges_sphere(les_faces_voisines(8),
864 les_faces_voisines(0),particle_id)*(valeurs_champ(les_faces_voisines(8)) -
865 valeurs_champ(les_faces_voisines(0))) / (xv(les_faces_voisines(8),1) -
866 xv(les_faces_voisines(0),1)) + compute_viscosity_edges_sphere(
867 les_faces_voisines(0),les_faces_voisines(10),particle_id)*(valeurs_champ(
868 les_faces_voisines(0)) - valeurs_champ(les_faces_voisines(10))) /
869 (xv(les_faces_voisines(0),1) - xv(les_faces_voisines(10),1)))
870 + 1./2.*( compute_viscosity_edges_sphere(les_faces_voisines(9),
871 les_faces_voisines(1),particle_id)*(valeurs_champ(les_faces_voisines(9)) -
872 valeurs_champ(les_faces_voisines(1))) / (xv(les_faces_voisines(9),1) -
873 xv(les_faces_voisines(1),1)) + compute_viscosity_edges_sphere(
874 les_faces_voisines(1),les_faces_voisines(11),particle_id)*(valeurs_champ(
875 les_faces_voisines(1)) - valeurs_champ(les_faces_voisines(11))) /
876 (xv(les_faces_voisines(1),1) - xv(les_faces_voisines(11),1))) );
877
878 gradU(elem,0,2) = 1./2.* ( 1./2.*( compute_viscosity_edges_sphere(les_faces_voisines(12),
879 les_faces_voisines(0),particle_id)*(valeurs_champ(les_faces_voisines(12)) -
880 valeurs_champ(les_faces_voisines(0))) / (xv(les_faces_voisines(12),2) -
881 xv(les_faces_voisines(0),2)) + compute_viscosity_edges_sphere(
882 les_faces_voisines(0),les_faces_voisines(13),particle_id)*(valeurs_champ(
883 les_faces_voisines(0)) - valeurs_champ(les_faces_voisines(13))) /
884 (xv(les_faces_voisines(0),2) - xv(les_faces_voisines(13),2)))
885 + 1./2.*( compute_viscosity_edges_sphere(les_faces_voisines(14),
886 les_faces_voisines(1),particle_id)*(valeurs_champ(les_faces_voisines(14)) -
887 valeurs_champ(les_faces_voisines(1))) / (xv(les_faces_voisines(14),2) -
888 xv(les_faces_voisines(1),2)) + compute_viscosity_edges_sphere(
889 les_faces_voisines(1),les_faces_voisines(15),particle_id)*(valeurs_champ(
890 les_faces_voisines(1)) - valeurs_champ(les_faces_voisines(15))) /
891 (xv(les_faces_voisines(1),2) - xv(les_faces_voisines(15),2))) );
892
893 gradU(elem,1,0) = 1./2.* ( 1./2.*( compute_viscosity_edges_sphere(les_faces_voisines(4),
894 les_faces_voisines(2),particle_id)*(valeurs_champ(les_faces_voisines(4)) -
895 valeurs_champ(les_faces_voisines(2))) / (xv(les_faces_voisines(4),0) -
896 xv(les_faces_voisines(2),0)) + compute_viscosity_edges_sphere(
897 les_faces_voisines(3),les_faces_voisines(5),particle_id)*(valeurs_champ(
898 les_faces_voisines(2)) - valeurs_champ(les_faces_voisines(5))) /
899 (xv(les_faces_voisines(2),0) - xv(les_faces_voisines(5),0)))
900 + 1./2.*( compute_viscosity_edges_sphere(les_faces_voisines(6),
901 les_faces_voisines(3),particle_id)*(valeurs_champ(les_faces_voisines(6)) -
902 valeurs_champ(les_faces_voisines(3))) / (xv(les_faces_voisines(6),0) -
903 xv(les_faces_voisines(3),0)) + compute_viscosity_edges_sphere(
904 les_faces_voisines(3),les_faces_voisines(7),particle_id)*(valeurs_champ(
905 les_faces_voisines(3)) - valeurs_champ(les_faces_voisines(7))) /
906 (xv(les_faces_voisines(3),0) - xv(les_faces_voisines(7),0))) );
907
908 gradU(elem,1,1) = ( (mu_p*mu_f/(mu_f-phase_indicator_function(elem_voisins(elem))*
909 (mu_f-mu_p)))*(valeurs_champ(les_faces_voisines(2)) - valeurs_champ(
910 les_faces_voisines(3))) / (xv(les_faces_voisines(2),1) - xv(les_faces_voisines(3),1)));
911
912 gradU(elem,1,2) = 1./2.* ( 1./2.*( compute_viscosity_edges_sphere(
913 les_faces_voisines(16),les_faces_voisines(2),particle_id)*(valeurs_champ(
914 les_faces_voisines(16)) - valeurs_champ(les_faces_voisines(2))) / (
915 xv(les_faces_voisines(16),2) - xv(les_faces_voisines(2),2)) +
916 compute_viscosity_edges_sphere(les_faces_voisines(2),les_faces_voisines(17),
917 particle_id)*(valeurs_champ(les_faces_voisines(2)) - valeurs_champ(
918 les_faces_voisines(17))) / (xv(les_faces_voisines(2),2) -
919 xv(les_faces_voisines(17),2)))
920 + 1./2.*( compute_viscosity_edges_sphere(les_faces_voisines(18),
921 les_faces_voisines(3),particle_id)*(valeurs_champ(les_faces_voisines(18))
922 - valeurs_champ(les_faces_voisines(3))) / (xv(les_faces_voisines(18),2) -
923 xv(les_faces_voisines(3),2)) + compute_viscosity_edges_sphere(
924 les_faces_voisines(3),les_faces_voisines(19),particle_id)*(valeurs_champ(
925 les_faces_voisines(3)) - valeurs_champ(les_faces_voisines(19))) /
926 (xv(les_faces_voisines(3),2) - xv(les_faces_voisines(19),2))));
927
928 gradU(elem,2,0) = 1./2.* ( 1./2.*( compute_viscosity_edges_sphere(les_faces_voisines(22),
929 les_faces_voisines(20),particle_id)*(valeurs_champ(les_faces_voisines(22)) -
930 valeurs_champ(les_faces_voisines(20))) / (xv(les_faces_voisines(22),0) -
931 xv(les_faces_voisines(20),0)) + compute_viscosity_edges_sphere(
932 les_faces_voisines(20),les_faces_voisines(24),particle_id)*(valeurs_champ(
933 les_faces_voisines(20)) - valeurs_champ(les_faces_voisines(24))) /
934 (xv(les_faces_voisines(20),0) - xv(les_faces_voisines(24),0)))
935 + 1./2.*( compute_viscosity_edges_sphere(les_faces_voisines(23),
936 les_faces_voisines(21),particle_id)*(valeurs_champ(les_faces_voisines(23)) -
937 valeurs_champ(les_faces_voisines(21))) / (xv(les_faces_voisines(23),0) -
938 xv(les_faces_voisines(21),0)) + compute_viscosity_edges_sphere(
939 les_faces_voisines(21),les_faces_voisines(25),particle_id)*(valeurs_champ(
940 les_faces_voisines(21)) - valeurs_champ(les_faces_voisines(25))) /
941 (xv(les_faces_voisines(21),0) - xv(les_faces_voisines(25),0))) );
942
943 gradU(elem,2,1) = 1./2.* ( 1./2.*( compute_viscosity_edges_sphere(les_faces_voisines(26),
944 les_faces_voisines(20),particle_id)*(valeurs_champ(les_faces_voisines(26)) -
945 valeurs_champ(les_faces_voisines(20))) / (xv(les_faces_voisines(26),1) -
946 xv(les_faces_voisines(20),1)) + compute_viscosity_edges_sphere(
947 les_faces_voisines(20),les_faces_voisines(28),particle_id)*(valeurs_champ(
948 les_faces_voisines(20)) - valeurs_champ(les_faces_voisines(28))) /
949 (xv(les_faces_voisines(20),1) - xv(les_faces_voisines(28),1)))
950 + 1./2.*( compute_viscosity_edges_sphere(les_faces_voisines(27),
951 les_faces_voisines(21),particle_id)*(valeurs_champ(les_faces_voisines(27)) -
952 valeurs_champ(les_faces_voisines(21))) / (xv(les_faces_voisines(27),1) -
953 xv(les_faces_voisines(21),1)) + compute_viscosity_edges_sphere(
954 les_faces_voisines(21),les_faces_voisines(29),particle_id)*(valeurs_champ(
955 les_faces_voisines(21)) - valeurs_champ(les_faces_voisines(29))) /
956 (xv(les_faces_voisines(21),1) - xv(les_faces_voisines(29),1))));
957
958 gradU(elem,2,2) = ( (mu_p*mu_f/(mu_f-phase_indicator_function(elem_voisins(elem))*
959 (mu_f-mu_p)))*(valeurs_champ(les_faces_voisines(20)) - valeurs_champ(
960 les_faces_voisines(21))) / (xv(les_faces_voisines(20),2) -
961 xv(les_faces_voisines(21),2)));
962 }
963
964 // trilinear interpolation of each component of the tensor
965 for (int i=0; i<dimension; i++)
966 {
967 for (int j=0; j<dimension; j++)
968 {
969 resu(fa7,i,j)=((1-zfact)*(
970 (1-yfact)*(
971 (1-xfact)*(gradU(0,i,j)) + xfact*(gradU(1,i,j))) +
972 yfact*((1-xfact)*(gradU(2,i,j)) + xfact*(gradU(3,i,j)))
973 ) +
974 zfact*(
975 (1-yfact)*((1-xfact)*(gradU(4,i,j)) + xfact*(gradU(5,i,j))) +
976 yfact*((1-xfact)*(gradU(6,i,j)) + xfact*(gradU(7,i,j))))
977 );
978 }
979 }
980 }
981 }
982 return 1;
983}
984
985int Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_gradU_elem(
986 const DoubleTab& valeurs_champ,
987 DoubleTab& coord,
988 DoubleTab& resu)
989{
990 Transport_Interfaces_FT_Disc& eq_transport = ptr_eq_transport_.valeur();
991 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
992 const Fluide_Diphasique& two_phase_fluid = eq_ns.fluide_diphasique();
993 const Maillage_FT_Disc& mesh = eq_transport.maillage_interface();
994 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, eq_ns.domaine_dis());
995 const DoubleTab& xv = domain_vdf.xv();
996 int nb_voisins=8;
997 const int id_fluid_phase=two_phase_fluid.get_id_fluid_phase();
998 const double mu_f =two_phase_fluid.fluide_phase(id_fluid_phase).
999 viscosite_dynamique().valeurs()(0, 0);
1000
1001 for (int fa7=0; fa7<coord.dimension(0); fa7++)
1002 {
1003 if (!mesh.facette_virtuelle(fa7))
1004 {
1005 // find the 8 neighboring elements where the stress tensor will be computed
1006 DoubleVect coord_elem_interp(dimension);
1007 for (int dim=0; dim<dimension; dim++) coord_elem_interp(dim)=coord(fa7,dim);
1008 IntVect elem_voisins(nb_voisins);
1009 find_neighboring_elements(coord_elem_interp,elem_voisins);
1010
1011 for (int i=0; i<nb_voisins; i++)
1012 {
1013 if (elem_voisins(i)<0)
1014 return 0;
1015 }
1016
1017 // compute distance between neighboring meshes in all direction to compute
1018 // interpolation coefficients
1019
1020 DoubleVect delta_i(dimension);
1021 delta_i(0) = fabs(domain_vdf.dist_elem(elem_voisins(0), elem_voisins(1), 0));
1022 delta_i(1) = fabs(domain_vdf.dist_elem(elem_voisins(0), elem_voisins(3), 1));
1023 delta_i(2) = fabs(domain_vdf.dist_elem(elem_voisins(0), elem_voisins(5), 2));
1024 DoubleVect coord_elem_0(dimension);
1025 for (int dim=0; dim<dimension; dim++) coord_elem_0(dim)=domain_vdf.xp(elem_voisins(0),dim);
1026
1027 double xfact=fabs((coord_elem_interp(0)-coord_elem_0(0))/delta_i(0));
1028 double yfact=fabs((coord_elem_interp(1)-coord_elem_0(1))/delta_i(1));
1029 double zfact=fabs((coord_elem_interp(2)-coord_elem_0(2))/delta_i(2));
1030
1031 // compute the velocity gradient tensor for each neighboring element
1032 DoubleTab gradU(nb_voisins, dimension, dimension);
1033 for (int elem=0; elem<nb_voisins; elem++)
1034 {
1035
1036 // SCHEMA EN 2D
1037 // AVANT ---- ---- --- --- --- --- --- --- --- --- AVANT
1038 // | | | | | | | | | |
1039 // 9 10
1040 // | | | | | | | | | |
1041 // ---- --- --- -5- -3- -6- --- --- ---
1042 // | | | | | | | | | |
1043 // GAUCHE 1 x 2 DROITE
1044 // | | | | | | | | | |
1045 // ---- --- --- -7- -4- -8- --- --- --- y
1046 // | | | | | | | | | | ^
1047 // 11 12 |
1048 // | | | | | | | | | | |
1049 // ARRIERE ---- --- --- --- --- --- --- --- --- ARRIERE o----->x
1050 // z
1051 // Evaluation du gradient sur les faces de normale x (choix arbitraire mais peut tuer une certaine symetrie)
1052 // Les termes + designent les faces juxtaposees selon +z
1053 // Les termes - designent les faces juxtaposees selon -z
1054 // Les termes 57 et 86 pour w designent les faces de normale z secantes aux faces 5,7 et 8,6
1055 // -- --
1056 // | u_2-u_1 1 u_9-u_11 u_10-u_12 1 u_1+ - u_1- u_2+ - u_2- |
1057 // | ---------- - (---------- + ----------- ) - ( ----------- + ----------- ) |
1058 // | x_2-x1 2 y_9-y_11 y_10-y_12 2 z_1+ - z_1- z_2+ - z_2- |
1059 // | |
1060 // grad U = | 1 v_5-v_6 v_7 - v_8 v_3-v_4 1 v_3+ -v_3- v_4+ - v_4- |
1061 // | - ( --------- + --------- ) ------- - ( ----------- + ----------- ) |
1062 // | 2 x_5-x_6 x_7 - x_8 y_3-y_4 2 z_3+ - z_3- z_4+ - z_4- |
1063 // | |
1064 // | 1 w57+ - w86+ w57- - w86- 1 w1112+ - w910+ w1112- -w910- w_34+ - w_34- |
1065 // | - ( ----------- + ----------- ) - ( -------------- + -------------- ) ------------- |
1066 // | 2 x57+ - x86+ x57- - x86- 2 y1112+ - y910+ y1112--y910- z_34+ - z_34- |
1067 // -- --
1068 //
1069 // Pour chaque face, il faut donc 30 faces voisines : 1-12, 1+, 1-, 2+, 2-, 3+, 3-, 4+, 4-, 34+, 34-, 57+, 57-, 86+, 86-, 1112+, 1112-, 910+, 910-
1070 // Pour plus de facilites, on numerote : 1+:13, 1-:14, 2+:15, 2-:16, 3+:17, 3-:18, 4+:19, 4-:20, 34+:21, 34-:22, 57+:23, 57-:24, 86+:25, 86-:26, 1112+:27, 1112-:28, 910+:29, 910-:30
1071
1072 // On identifie les faces voisines pour le calcul des composantes du tenseur
1073 int nb_faces_voisines=30;
1074 IntVect les_faces_voisines(nb_faces_voisines);
1075 int elem_=elem_voisins(elem);
1076 int elem_gauche = face_voisins_for_interp(elem_faces_for_interp(elem_,0),0);
1077 int elem_droite = face_voisins_for_interp(elem_faces_for_interp(elem_,0+dimension),1);
1078 int elem_haut=face_voisins_for_interp(elem_faces_for_interp(elem_,2+dimension), 1);
1079 int elem_bas=face_voisins_for_interp(elem_faces_for_interp(elem_,2), 0);
1080 int elem_avant=face_voisins_for_interp(elem_faces_for_interp(elem_,1+dimension),1);
1081 int elem_arriere=face_voisins_for_interp(elem_faces_for_interp(elem_,1),0);
1082
1083 les_faces_voisines(0)=elem_faces_for_interp(elem_,0);
1084 les_faces_voisines(1)=elem_faces_for_interp(elem_,0+dimension);
1085 les_faces_voisines(2)=elem_faces_for_interp(elem_,1);
1086 les_faces_voisines(3)=elem_faces_for_interp(elem_,1+dimension);
1087 les_faces_voisines(4)=elem_faces_for_interp(elem_gauche,1+dimension);
1088 les_faces_voisines(5)=elem_faces_for_interp(elem_droite,1+dimension);
1089 les_faces_voisines(6)=elem_faces_for_interp(elem_gauche,1);
1090 les_faces_voisines(7)=elem_faces_for_interp(elem_droite,1);
1091 les_faces_voisines(8)=elem_faces_for_interp(elem_avant,0);
1092 les_faces_voisines(9)=elem_faces_for_interp(elem_avant,0+dimension);
1093 les_faces_voisines(10)=elem_faces_for_interp(elem_arriere,0);
1094 les_faces_voisines(11)=elem_faces_for_interp(elem_arriere,0+dimension);
1095 les_faces_voisines(12)=elem_faces_for_interp(elem_haut,0);
1096 les_faces_voisines(13)=elem_faces_for_interp(elem_bas,0);
1097 les_faces_voisines(14)=elem_faces_for_interp(elem_haut,0+dimension);
1098 les_faces_voisines(15)=elem_faces_for_interp(elem_bas,0+dimension);
1099 les_faces_voisines(16)=elem_faces_for_interp(elem_haut,1+dimension);
1100 les_faces_voisines(17)=elem_faces_for_interp(elem_bas,1+dimension);
1101 les_faces_voisines(18)=elem_faces_for_interp(elem_haut,1);
1102 les_faces_voisines(19)=elem_faces_for_interp(elem_bas,1);
1103 les_faces_voisines(20)=elem_faces_for_interp(elem_,2);
1104 les_faces_voisines(21)=elem_faces_for_interp(elem_,2+dimension);
1105 les_faces_voisines(22)=elem_faces_for_interp(elem_gauche,2);
1106 les_faces_voisines(23)=elem_faces_for_interp(elem_gauche,2+dimension);
1107 les_faces_voisines(24)=elem_faces_for_interp(elem_droite,2);
1108 les_faces_voisines(25)=elem_faces_for_interp(elem_droite,2+dimension);
1109 les_faces_voisines(26)=elem_faces_for_interp(elem_arriere,2);
1110 les_faces_voisines(27)=elem_faces_for_interp(elem_arriere,2+dimension);
1111 les_faces_voisines(28)=elem_faces_for_interp(elem_avant,2);
1112 les_faces_voisines(29)=elem_faces_for_interp(elem_avant,2+dimension);
1113
1114 for (int i=0; i<nb_faces_voisines; i++)
1115 {
1116 if (les_faces_voisines(i)<0)
1117 return 0;
1118 }
1119
1120 gradU(elem,0,0) = ( (valeurs_champ(les_faces_voisines(1)) - valeurs_champ(
1121 les_faces_voisines(0))) / (xv(les_faces_voisines(1),0) - xv(les_faces_voisines(0),0)));
1122
1123 gradU(elem,0,1) = 1./2.*( (valeurs_champ(les_faces_voisines(8)) - valeurs_champ(
1124 les_faces_voisines(10))) / (xv(les_faces_voisines(8),1) -
1125 xv(les_faces_voisines(10),1)) + (valeurs_champ(les_faces_voisines(9)) -
1126 valeurs_champ(les_faces_voisines(11))) / (xv(les_faces_voisines(9),1) -
1127 xv(les_faces_voisines(11),1)));
1128
1129 gradU(elem,0,2) = 1./2.*( (valeurs_champ(les_faces_voisines(12)) - valeurs_champ(
1130 les_faces_voisines(13))) / (xv(les_faces_voisines(12),2) -
1131 xv(les_faces_voisines(13),2)) + (valeurs_champ(les_faces_voisines(14)) -
1132 valeurs_champ(les_faces_voisines(15))) / (xv(les_faces_voisines(14),2) -
1133 xv(les_faces_voisines(15),2)));
1134
1135 gradU(elem,1,0) = 1./2.*( (valeurs_champ(les_faces_voisines(4)) - valeurs_champ(
1136 les_faces_voisines(5))) / (xv(les_faces_voisines(4),0) -
1137 xv(les_faces_voisines(5),0)) + (valeurs_champ(les_faces_voisines(6)) -
1138 valeurs_champ(les_faces_voisines(7))) / (xv(les_faces_voisines(6),0) -
1139 xv(les_faces_voisines(7),0)));
1140
1141 gradU(elem,1,1) = ( (valeurs_champ(les_faces_voisines(2)) - valeurs_champ(
1142 les_faces_voisines(3))) / (xv(les_faces_voisines(2),1) -
1143 xv(les_faces_voisines(3),1)));
1144
1145 gradU(elem,1,2) = 1./2.*( (valeurs_champ(les_faces_voisines(16)) - valeurs_champ(
1146 les_faces_voisines(17))) / (xv(les_faces_voisines(16),2) -
1147 xv(les_faces_voisines(17),2)) + (valeurs_champ(les_faces_voisines(18)) -
1148 valeurs_champ(les_faces_voisines(19))) / (xv(les_faces_voisines(18),2) -
1149 xv(les_faces_voisines(19),2)));
1150
1151 gradU(elem,2,0) = 1./2.*( (valeurs_champ(les_faces_voisines(22)) - valeurs_champ(
1152 les_faces_voisines(24))) / (xv(les_faces_voisines(22),0) -
1153 xv(les_faces_voisines(24),0)) + (valeurs_champ(les_faces_voisines(23)) -
1154 valeurs_champ(les_faces_voisines(25))) / (xv(les_faces_voisines(23),0) -
1155 xv(les_faces_voisines(25),0)));
1156
1157 gradU(elem,2,1) = 1./2.*( (valeurs_champ(les_faces_voisines(26)) - valeurs_champ(
1158 les_faces_voisines(28))) / (xv(les_faces_voisines(26),1) -
1159 xv(les_faces_voisines(28),1)) + (valeurs_champ(les_faces_voisines(27)) -
1160 valeurs_champ(les_faces_voisines(29))) / (xv(les_faces_voisines(27),1) -
1161 xv(les_faces_voisines(29),1)));
1162
1163 gradU(elem,2,2) = ( (valeurs_champ(les_faces_voisines(20)) - valeurs_champ(
1164 les_faces_voisines(21))) / (xv(les_faces_voisines(20),2) -
1165 xv(les_faces_voisines(21),2)));
1166 }
1167
1168 // On fait une interpolation trilineaire de chaque composante du tenseur
1169 for (int i=0; i<dimension; i++)
1170 {
1171 for (int j=0; j<dimension; j++)
1172 {
1173 resu(fa7,i,j)=mu_f*((1-zfact)*(
1174 (1-yfact)*(
1175 (1-xfact)*(gradU(0,i,j)) + xfact*(gradU(1,i,j))) +
1176 yfact*((1-xfact)*(gradU(2,i,j)) + xfact*(gradU(3,i,j)))
1177 ) +
1178 zfact*(
1179 (1-yfact)*((1-xfact)*(gradU(4,i,j)) + xfact*(gradU(5,i,j))) +
1180 yfact*((1-xfact)*(gradU(6,i,j)) + xfact*(gradU(7,i,j))))
1181 );
1182 }
1183 }
1184 }
1185 }
1186 return 1;
1187}
1188
1189int Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_face_sommets(
1190 const DoubleTab& valeurs_champ,
1191 DoubleTab& coord,
1192 DoubleTab& resu)
1193{
1194 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
1195 const Transport_Interfaces_FT_Disc& eq_transport = ptr_eq_transport_.valeur();
1196 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, eq_ns.domaine_dis());
1197 const Maillage_FT_Disc& mesh = eq_transport.maillage_interface();
1198 const ArrOfInt& elem = mesh.sommet_elem();
1199 const DoubleTab& pos = mesh.sommets();
1200 const int nb_pos_tot = pos.dimension(0);
1201 int nb_voisins=8;
1202 for (int som=0; som<nb_pos_tot; som++)
1203 {
1204 const int element = elem[som];
1205 if (element >= 0) // real vertice ?
1206 {
1207 // find the neighboring faces that will be used for interpolation
1208 DoubleVect coord_elem_interp(dimension); // coordinate of the fa7 gravity center
1209 for (int dim=0; dim<dimension; dim++)
1210 coord_elem_interp(dim)=coord(som,dim);
1211 IntTab faces_voisines(dimension,nb_voisins);
1212 find_neighboring_faces_xyz(coord_elem_interp,faces_voisines);
1213 for (int dim=0; dim<dimension; dim++)
1214 {
1215 for (int i=0; i<nb_voisins; i++)
1216 {
1217 if (faces_voisines(dim,i)<0)
1218 return 0; // if one face could not be computed (access problem : joint zone, wall),
1219 // we return 0 and the fa7 is not considered in the computation
1220 }
1221 }
1222
1223 DoubleVect xfact(dimension);
1224 DoubleVect yfact(dimension);
1225 DoubleVect zfact(dimension);
1226 DoubleTab Delta_(dimension);
1227 Delta_(0)=fabs(domain_vdf.dist_face(faces_voisines(0,0),faces_voisines(0,1),0));
1228 Delta_(1)=fabs(domain_vdf.dist_face(faces_voisines(0,0),faces_voisines(0,2),1));
1229 Delta_(2)=fabs(domain_vdf.dist_face(faces_voisines(0,0),faces_voisines(0,4),2));
1230
1231 for (int dim=0; dim<dimension; dim++)
1232 {
1233 xfact(dim)=fabs((coord(som,0)-domain_vdf.xv(faces_voisines(dim,0),0))/Delta_(0));
1234 yfact(dim)=fabs((coord(som,1)-domain_vdf.xv(faces_voisines(dim,0),1))/Delta_(1));
1235 zfact(dim)=fabs((coord(som,2)-domain_vdf.xv(faces_voisines(dim,0),2))/Delta_(2));
1236 }
1237
1238 for (int dim=0; dim<dimension; dim++)
1239 {
1240 if (xfact(dim)>1) Cerr << "xfact > 1 pour le point " <<coord(som,0)<< " " << coord(som,1)<< " " << coord(som,2) <<finl;
1241 if (yfact(dim)>1) Cerr << "yfact > 1 pour le point " <<coord(som,0)<< " " << coord(som,1)<< " " << coord(som,2) <<finl;
1242 if (zfact(dim)>1) Cerr << "zfact > 1 pour le point " <<coord(som,0)<< " " << coord(som,1)<< " " << coord(som,2) <<finl;
1243 }
1244 for (int dim=0; dim<dimension; dim++)
1245 {
1246 resu(som,dim)=(1.-zfact(dim))*(
1247 (1.-yfact(dim))*(
1248 (1.-xfact(dim))*(valeurs_champ(faces_voisines(dim,0))) +
1249 xfact(dim)*(valeurs_champ(faces_voisines(dim,1)))
1250 ) +
1251 yfact(dim)*(
1252 (1.-xfact(dim))*(valeurs_champ(faces_voisines(dim,2))) +
1253 xfact(dim)*(valeurs_champ(faces_voisines(dim,3))))
1254 ) +
1255 zfact(dim)*(
1256 (1.-yfact(dim))*(
1257 (1.-xfact(dim))*(valeurs_champ(faces_voisines(dim,4))) +
1258 xfact(dim)*(valeurs_champ(faces_voisines(dim,5)))
1259 ) +
1260 yfact(dim)*(
1261 (1.-xfact(dim))*(valeurs_champ(faces_voisines(dim,6))) +
1262 xfact(dim)*(valeurs_champ(faces_voisines(dim,7))))
1263 );
1264 }
1265 }
1266 }
1267
1268 return 1;
1269}
1270
1271double Post_Processing_Hydrodynamic_Forces::find_neighboring_elements(
1272 DoubleVect& coord_elem_interp,
1273 IntVect& neighboring_elements,
1274 const int sauv_list_P1,
1275 const int num_fa7)
1276{
1277 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
1278 Transport_Interfaces_FT_Disc& eq_transport = ptr_eq_transport_.valeur();
1279 const DoubleTab& phase_indicator_function = eq_transport.get_indicatrice().valeurs();
1280 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, eq_ns.domaine_dis());
1281 const Domaine& domain = domain_vdf.domaine();
1282
1283 DoubleVect coord_elem_eulerian(dimension);
1284 int elem_eulerien=domain.chercher_elements(coord_elem_interp(0),
1285 coord_elem_interp(1),
1286 coord_elem_interp(2));
1287 if (sauv_list_P1)
1288 list_elem_P1_(num_fa7,0)=elem_eulerien;
1289 if (elem_eulerien<0)
1290 {
1291 neighboring_elements=-1;
1292 return -1;
1293 }
1294
1295 for (int dim=0; dim<dimension; dim++)
1296 coord_elem_eulerian(dim)=domain_vdf.xp(elem_eulerien,dim);
1297 IntVect direction_interp(dimension); // Pour chaque direction, on regarde de quel cote de l'element le point se trouve. 0 : a gauche (pos_i_point<=pos_i_elem), 1 : a droite(pos_i_point>pos_i_elem)
1298 IntVect faces_elem_interp(2*dimension);
1299 for (int dim=0; dim<dimension; dim++)
1300 {
1301 faces_elem_interp(dim)=elem_faces_for_interp(elem_eulerien,dim);
1302 faces_elem_interp(dimension+dim)=elem_faces_for_interp(elem_eulerien,dimension+dim);
1303
1304 if (coord_elem_interp(dim)<=coord_elem_eulerian(dim))
1305 direction_interp(dim)=0;
1306 else
1307 direction_interp(dim)=1;
1308 }
1309
1310 // The 8 elements form a large cube with 2 cubes in each direction (1 cube = 1 elem).
1311 // Let's consider a cube with side 1. Each vertice represent a neighboring element.
1312 // We fill the neighboring list beginning with coordinates
1313 // z=0, y=0 then z=0, y=1, then z=1, y=0 then z=1, y=1
1314 // Thus, we have (0,0,0) ; (1,0,0) ; (0,1,0) ; (1,1,0) ; (0,0,+-1) ; (1,0,+-1) ; (0,1,+-1) ; (1,1,+-1)
1315 if (direction_interp(0)==1)
1316 {
1317 if (direction_interp(1)==1)
1318 {
1319 neighboring_elements(0)=elem_eulerien;
1320 neighboring_elements(1)=face_voisins_for_interp(faces_elem_interp(0+dimension),1);
1321 neighboring_elements(2)=face_voisins_for_interp(faces_elem_interp(1+dimension),1);
1322 neighboring_elements(3)=face_voisins_for_interp(elem_faces_for_interp(
1323 neighboring_elements(1),1+dimension),1);
1324 }
1325 else
1326 {
1327 neighboring_elements(0)=face_voisins_for_interp(faces_elem_interp(1),0);
1328 neighboring_elements(1)=face_voisins_for_interp(elem_faces_for_interp(
1329 neighboring_elements(0),0+dimension),1);
1330 neighboring_elements(2)=elem_eulerien;
1331 neighboring_elements(3)=face_voisins_for_interp(faces_elem_interp(0+dimension),1);
1332 }
1333 }
1334 else
1335 {
1336 if (direction_interp(1)==1)
1337 {
1338 neighboring_elements(0)=face_voisins_for_interp(faces_elem_interp(0),0);
1339 neighboring_elements(1)=elem_eulerien;
1340 neighboring_elements(2)=face_voisins_for_interp(elem_faces_for_interp(
1341 neighboring_elements(0),1+dimension),1);
1342 neighboring_elements(3)=face_voisins_for_interp(faces_elem_interp(1+dimension),1);
1343 }
1344 else
1345 {
1346 neighboring_elements(1)=face_voisins_for_interp(faces_elem_interp(1),0);
1347 neighboring_elements(2)=face_voisins_for_interp(faces_elem_interp(0),0);
1348 neighboring_elements(3)=elem_eulerien;
1349 neighboring_elements(0)=face_voisins_for_interp(elem_faces_for_interp(
1350 neighboring_elements(1),0),0);
1351 }
1352 }
1353
1354 if (direction_interp(2)==1)
1355 {
1356 neighboring_elements(4)=face_voisins_for_interp(elem_faces_for_interp(
1357 neighboring_elements(0),2+dimension),1);
1358 neighboring_elements(5)=face_voisins_for_interp(elem_faces_for_interp(
1359 neighboring_elements(1),2+dimension),1);
1360 neighboring_elements(6)=face_voisins_for_interp(elem_faces_for_interp(
1361 neighboring_elements(2),2+dimension),1);
1362 neighboring_elements(7)=face_voisins_for_interp(elem_faces_for_interp(
1363 neighboring_elements(3),2+dimension),1);
1364 }
1365 else
1366 {
1367 neighboring_elements(4)=face_voisins_for_interp(elem_faces_for_interp(
1368 neighboring_elements(0),2),0);
1369 neighboring_elements(5)=face_voisins_for_interp(elem_faces_for_interp(
1370 neighboring_elements(1),2),0);
1371 neighboring_elements(6)=face_voisins_for_interp(elem_faces_for_interp(
1372 neighboring_elements(2),2),0);
1373 neighboring_elements(7)=face_voisins_for_interp(elem_faces_for_interp(
1374 neighboring_elements(3),2),0);
1375
1376 for (int i=0; i<4; i++)
1377 {
1378 int tmp=neighboring_elements(i);
1379 neighboring_elements(i)=neighboring_elements(i+4);
1380 neighboring_elements(i+4)=tmp;
1381 }
1382 }
1383 return (phase_indicator_function(elem_eulerien));
1384}
1385
1386void Post_Processing_Hydrodynamic_Forces::find_neighboring_faces(
1387 DoubleVect& coord_elem_interp,
1388 IntVect& neighboring_faces,
1389 int orientation) const
1390{
1391 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
1392 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, eq_ns.domaine_dis());
1393
1394 DoubleVect coord_elem_eulerien(dimension);
1395
1396 int elem_eulerien=domain_vdf.domaine().chercher_elements(coord_elem_interp(0),
1397 coord_elem_interp(1),
1398 coord_elem_interp(2));
1399 if (elem_eulerien<0)
1400 {
1401 neighboring_faces=-1;
1402 return;
1403 }
1404 for (int dim=0; dim<dimension; dim++)
1405 coord_elem_eulerien(dim)=domain_vdf.xp(elem_eulerien,dim);
1406 // In each direction, we look at which side of the element the point is.
1407 // 0 : left (pos_i_point<=pos_i_elem), 1 : right (pos_i_point>pos_i_elem)
1408 IntVect direction_interp(dimension);
1409 IntVect faces_elem_interp(2*dimension);
1410
1411 for (int dim=0; dim<dimension; dim++)
1412 {
1413 faces_elem_interp(dim)=elem_faces_for_interp(elem_eulerien,dim);
1414 faces_elem_interp(dimension+dim)=elem_faces_for_interp(elem_eulerien,dimension+dim);
1415
1416 if (coord_elem_interp(dim)<=coord_elem_eulerien(dim))
1417 direction_interp(dim)=0;
1418 else
1419 direction_interp(dim)=1;
1420 }
1421
1422 if (orientation==0)
1423 {
1424 if (direction_interp(1)==1)
1425 {
1426 neighboring_faces(0)=faces_elem_interp(orientation);
1427 neighboring_faces(1)=faces_elem_interp(orientation+dimension);
1428 neighboring_faces(2)=elem_faces_for_interp(face_voisins_for_interp(
1429 faces_elem_interp(1+dimension),1),orientation);
1430 neighboring_faces(3)=elem_faces_for_interp(face_voisins_for_interp(
1431 faces_elem_interp(1+dimension),1),orientation+dimension);
1432 }
1433 else
1434 {
1435 neighboring_faces(0)=elem_faces_for_interp(face_voisins_for_interp(
1436 faces_elem_interp(1),0),orientation);
1437 neighboring_faces(1)=elem_faces_for_interp(face_voisins_for_interp(
1438 faces_elem_interp(1),0),orientation+dimension);
1439 neighboring_faces(2)=faces_elem_interp(orientation);
1440 neighboring_faces(3)=faces_elem_interp(orientation+dimension);
1441 }
1442 }
1443 if (orientation==1)
1444 {
1445 if (direction_interp(0)==1)
1446 {
1447 neighboring_faces(0)=faces_elem_interp(orientation);
1448 neighboring_faces(1)=elem_faces_for_interp(face_voisins_for_interp(
1449 faces_elem_interp(0+dimension),1),orientation);
1450 neighboring_faces(2)=faces_elem_interp(orientation+dimension);
1451 neighboring_faces(3)=elem_faces_for_interp(face_voisins_for_interp(
1452 faces_elem_interp(0+dimension),1),orientation+dimension);
1453 }
1454 else
1455 {
1456 neighboring_faces(0)=elem_faces_for_interp(face_voisins_for_interp(
1457 faces_elem_interp(0),0),orientation);
1458 neighboring_faces(1)=faces_elem_interp(orientation);
1459 neighboring_faces(2)=elem_faces_for_interp(face_voisins_for_interp(
1460 faces_elem_interp(0),0),orientation+dimension);
1461 neighboring_faces(3)=faces_elem_interp(orientation+dimension);
1462 }
1463 }
1464 if (direction_interp(2)==1 && orientation!=2)
1465 {
1466 neighboring_faces(4)=elem_faces_for_interp(face_voisins_for_interp(
1467 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(0),1),
1468 2+dimension),1),orientation);
1469 neighboring_faces(5)=elem_faces_for_interp(face_voisins_for_interp(
1470 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(1),1),
1471 2+dimension),1),orientation);
1472 neighboring_faces(6)=elem_faces_for_interp(face_voisins_for_interp(
1473 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(2),1),
1474 2+dimension),1),orientation);
1475 neighboring_faces(7)=elem_faces_for_interp(face_voisins_for_interp(
1476 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(3),1),
1477 2+dimension),1),orientation);
1478 }
1479 else if (direction_interp(2)==0 && orientation!=2)
1480 {
1481 neighboring_faces(4)=elem_faces_for_interp(face_voisins_for_interp(
1482 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(0),1),
1483 2),0),orientation);
1484 neighboring_faces(5)=elem_faces_for_interp(face_voisins_for_interp(
1485 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(1),1),
1486 2),0),orientation);
1487 neighboring_faces(6)=elem_faces_for_interp(face_voisins_for_interp(
1488 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(2),1),
1489 2),0),orientation);
1490 neighboring_faces(7)=elem_faces_for_interp(face_voisins_for_interp(
1491 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(3),1)
1492 ,2),0),orientation);
1493
1494 for (int i=0; i<4; i++)
1495 {
1496 int tmp=neighboring_faces(i);
1497 neighboring_faces(i)=neighboring_faces(i+4);
1498 neighboring_faces(i+4)=tmp;
1499 }
1500 }
1501
1502 if (orientation==2)
1503 {
1504 if (direction_interp(0)==1)
1505 {
1506 neighboring_faces(0)=faces_elem_interp(orientation);
1507 neighboring_faces(1)=elem_faces_for_interp(face_voisins_for_interp(
1508 faces_elem_interp(0+dimension),1),orientation);
1509 neighboring_faces(4)=faces_elem_interp(orientation+dimension);
1510 neighboring_faces(5)=elem_faces_for_interp(face_voisins_for_interp(
1511 faces_elem_interp(0+dimension),1),orientation+dimension);
1512 }
1513 else
1514 {
1515 neighboring_faces(0)=elem_faces_for_interp(face_voisins_for_interp(
1516 faces_elem_interp(0),0),orientation);
1517 neighboring_faces(1)=faces_elem_interp(orientation);
1518 neighboring_faces(4)=elem_faces_for_interp(face_voisins_for_interp(
1519 faces_elem_interp(0),0),orientation+dimension);
1520 neighboring_faces(5)=faces_elem_interp(orientation+dimension);
1521
1522 }
1523
1524 if (direction_interp(1)==1)
1525 {
1526 neighboring_faces(2)=elem_faces_for_interp(face_voisins_for_interp(
1527 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(0),1),
1528 1+dimension),1),orientation);
1529 neighboring_faces(3)=elem_faces_for_interp(face_voisins_for_interp(
1530 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(1),1),
1531 1+dimension),1),orientation);
1532 neighboring_faces(6)=elem_faces_for_interp(face_voisins_for_interp(
1533 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(4),1),
1534 1+dimension),1),orientation);
1535 neighboring_faces(7)=elem_faces_for_interp(face_voisins_for_interp(
1536 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(5),1),
1537 1+dimension),1),orientation);
1538 }
1539 else
1540 {
1541 neighboring_faces(2)=elem_faces_for_interp(face_voisins_for_interp(
1542 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(0),1),
1543 1),0),orientation);
1544 neighboring_faces(3)=elem_faces_for_interp(face_voisins_for_interp(
1545 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(1),1),
1546 1),0),orientation);
1547 neighboring_faces(6)=elem_faces_for_interp(face_voisins_for_interp(
1548 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(4),1),
1549 1),0),orientation);
1550 neighboring_faces(7)=elem_faces_for_interp(face_voisins_for_interp(
1551 elem_faces_for_interp(face_voisins_for_interp(neighboring_faces(5),1),
1552 1),0),orientation);
1553
1554 int tmp0=neighboring_faces(0);
1555 int tmp1=neighboring_faces(1);
1556 int tmp4=neighboring_faces(4);
1557 int tmp5=neighboring_faces(5);
1558
1559 neighboring_faces(0)=neighboring_faces(2);
1560 neighboring_faces(1)=neighboring_faces(3);
1561 neighboring_faces(4)=neighboring_faces(6);
1562 neighboring_faces(5)=neighboring_faces(7);
1563
1564 neighboring_faces(2)=tmp0;
1565 neighboring_faces(3)=tmp1;
1566 neighboring_faces(6)=tmp4;
1567 neighboring_faces(7)=tmp5;
1568 }
1569 }
1570}
1571
1572void Post_Processing_Hydrodynamic_Forces::find_neighboring_faces_xyz(
1573 DoubleVect& coord_elem_interp,
1574 IntTab& neighboring_faces) const
1575{
1576 int nb_neighboring_faces=8;
1577 IntVect neighboring_faces_x(nb_neighboring_faces);
1578 IntVect neighboring_faces_y(nb_neighboring_faces);
1579 IntVect neighboring_faces_z(nb_neighboring_faces);
1580
1581 find_neighboring_faces(coord_elem_interp,neighboring_faces_x,0);
1582 find_neighboring_faces(coord_elem_interp,neighboring_faces_y,1);
1583 find_neighboring_faces(coord_elem_interp,neighboring_faces_z,2);
1584
1585 for (int i=0; i<8; i++)
1586 {
1587 neighboring_faces(0,i)=neighboring_faces_x(i);
1588 neighboring_faces(1,i)=neighboring_faces_y(i);
1589 neighboring_faces(2,i)=neighboring_faces_z(i);
1590 }
1591}
1592
1593// see Vincent et al., 2014
1594double Post_Processing_Hydrodynamic_Forces::compute_viscosity_edges_sphere(int face1, int face2, int particle_id)
1595{
1596 Transport_Interfaces_FT_Disc& eq_transport = ptr_eq_transport_.valeur();
1597 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
1598 const Fluide_Diphasique& two_phase_fluid = eq_ns.fluide_diphasique();
1599
1600 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, eq_ns.domaine_dis());
1601 const DoubleTab& xp = domain_vdf.xp();
1602 const DoubleTab& particles_position=eq_transport.get_particles_position();
1603
1604 const int id_fluid_phase=two_phase_fluid.get_id_fluid_phase();
1605 const Solid_Particle_base& solid_particle=ref_cast(Solid_Particle_base,two_phase_fluid.fluide_phase(id_fluid_phase));
1606 const double& effective_radius=solid_particle.get_equivalent_radius();
1607 const double mu_f =two_phase_fluid.fluide_phase(id_fluid_phase).
1608 viscosite_dynamique().valeurs()(0, 0);
1609 const double mu_p = two_phase_fluid.fluide_phase(1-id_fluid_phase).
1610 viscosite_dynamique().valeurs()(0, 0);
1611
1612 int elem1 = face_voisins_for_interp(face1,0);
1613 int elem2 = face_voisins_for_interp(face1,1);
1614 int elem3 = face_voisins_for_interp(face2,0);
1615 int elem4 = face_voisins_for_interp(face2,1);
1616
1617 double indicatrice_arete=0;
1618 int pvx=10;
1619 int pvy=10;
1620 int pvz=10;
1621 double x_min,y_min,z_min,x_max;
1622 double x_cg=particles_position(particle_id,0);
1623 double y_cg=particles_position(particle_id,1);
1624 double z_cg=particles_position(particle_id,2);
1625
1626 x_min=std::min(xp(elem1,0),xp(elem2,0));
1627 x_min=std::min(x_min,xp(elem3,0));
1628 x_min=std::min(x_min,xp(elem4,0));
1629 x_max=std::max(xp(elem1,0),xp(elem2,0));
1630 x_max=std::max(x_min,xp(elem3,0));
1631 x_max=std::max(x_min,xp(elem4,0));
1632 y_min=std::min(xp(elem1,1),xp(elem2,1));
1633 y_min=std::min(y_min,xp(elem3,1));
1634 y_min=std::min(y_min,xp(elem4,1));
1635 z_min=std::min(xp(elem1,2),xp(elem2,2));
1636 z_min=std::min(z_min,xp(elem3,2));
1637 z_min=std::min(z_min,xp(elem4,2));
1638
1639 double delta_x=x_max-x_min;
1640
1641 for (int i=0; i<pvx; i++)
1642 {
1643 for (int j=0; j<pvy; j++)
1644 {
1645 for (int k=0; k<pvz; k++)
1646 {
1647 double x=x_min+i*delta_x/pvx;
1648 double y=y_min+i*delta_x/pvy;
1649 double z=z_min+i*delta_x/pvz;
1650 if (sqrt(pow(x-x_cg,2)+pow(y-y_cg,2)+pow(z-z_cg,2))>effective_radius)
1651 indicatrice_arete+=1;
1652 }
1653 }
1654 }
1655 indicatrice_arete/=(pvx*pvy*pvz);
1656
1657 return(mu_p*mu_f/(mu_f-indicatrice_arete*(mu_f-mu_p)));
1658}
1659
1660void Post_Processing_Hydrodynamic_Forces::resize_sigma(int nb_fa7)
1661{
1662 sigma_xx_fa7_.resize(nb_fa7);
1663 sigma_xy_fa7_.resize(nb_fa7);
1664 sigma_xz_fa7_.resize(nb_fa7);
1665 sigma_yx_fa7_.resize(nb_fa7);
1666 sigma_yy_fa7_.resize(nb_fa7);
1667 sigma_yz_fa7_.resize(nb_fa7);
1668 sigma_zx_fa7_.resize(nb_fa7);
1669 sigma_zy_fa7_.resize(nb_fa7);
1670 sigma_zz_fa7_.resize(nb_fa7);
1671}
1672
1673void Post_Processing_Hydrodynamic_Forces::resize_gradU_P1(int nb_fa7)
1674{
1675 dUdx_P1_.resize(nb_fa7);
1676 dUdy_P1_.resize(nb_fa7);
1677 dUdz_P1_.resize(nb_fa7);
1678 dVdx_P1_.resize(nb_fa7);
1679 dVdy_P1_.resize(nb_fa7);
1680 dVdz_P1_.resize(nb_fa7);
1681 dWdx_P1_.resize(nb_fa7);
1682 dWdy_P1_.resize(nb_fa7);
1683 dWdz_P1_.resize(nb_fa7);
1684}
1685
1686void Post_Processing_Hydrodynamic_Forces::resize_gradU_P2(int nb_fa7)
1687{
1688 dUdx_P2_.resize(nb_fa7);
1689 dUdy_P2_.resize(nb_fa7);
1690 dUdz_P2_.resize(nb_fa7);
1691 dVdx_P2_.resize(nb_fa7);
1692 dVdy_P2_.resize(nb_fa7);
1693 dVdz_P2_.resize(nb_fa7);
1694 dWdx_P2_.resize(nb_fa7);
1695 dWdy_P2_.resize(nb_fa7);
1696 dWdz_P2_.resize(nb_fa7);
1697}
1698
1699void Post_Processing_Hydrodynamic_Forces::fill_gradU_P1(int fa7, DoubleTab gradU_P1)
1700{
1701 dUdx_P1_(fa7)=gradU_P1(fa7,0,0);
1702 dUdy_P1_(fa7)=gradU_P1(fa7,0,1);
1703 dUdz_P1_(fa7)=gradU_P1(fa7,0,2);
1704 dVdx_P1_(fa7)=gradU_P1(fa7,1,0);
1705 dVdy_P1_(fa7)=gradU_P1(fa7,1,1);
1706 dVdz_P1_(fa7)=gradU_P1(fa7,1,2);
1707 dWdx_P1_(fa7)=gradU_P1(fa7,2,0);
1708 dWdy_P1_(fa7)=gradU_P1(fa7,2,1);
1709 dWdz_P1_(fa7)=gradU_P1(fa7,2,2);
1710}
1711
1712void Post_Processing_Hydrodynamic_Forces::fill_gradU_P2(int fa7, DoubleTab gradU_P2)
1713{
1714 dUdx_P2_(fa7)=gradU_P2(fa7,0,0);
1715 dUdy_P2_(fa7)=gradU_P2(fa7,0,1);
1716 dUdz_P2_(fa7)=gradU_P2(fa7,0,2);
1717 dVdx_P2_(fa7)=gradU_P2(fa7,1,0);
1718 dVdy_P2_(fa7)=gradU_P2(fa7,1,1);
1719 dVdz_P2_(fa7)=gradU_P2(fa7,1,2);
1720 dWdx_P2_(fa7)=gradU_P2(fa7,2,0);
1721 dWdy_P2_(fa7)=gradU_P2(fa7,2,1);
1722 dWdz_P2_(fa7)=gradU_P2(fa7,2,2);
1723}
1724
1725void Post_Processing_Hydrodynamic_Forces::fill_sigma(int fa7, Matrice_Dense stress_tensor)
1726{
1727 sigma_xx_fa7_(fa7)=stress_tensor(0,0);
1728 sigma_xy_fa7_(fa7)=stress_tensor(0,1);
1729 sigma_xz_fa7_(fa7)=stress_tensor(0,2);
1730 sigma_yx_fa7_(fa7)=stress_tensor(1,0);
1731 sigma_yy_fa7_(fa7)=stress_tensor(1,1);
1732 sigma_yz_fa7_(fa7)=stress_tensor(1,2);
1733 sigma_zx_fa7_(fa7)=stress_tensor(2,0);
1734 sigma_zy_fa7_(fa7)=stress_tensor(2,1);
1735 sigma_zz_fa7_(fa7)=stress_tensor(2,2);
1736}
1737
1738void Post_Processing_Hydrodynamic_Forces::resize_and_init_tables(int nb_particles_tot)
1739{
1740 total_pressure_force_.resize(nb_particles_tot,dimension);
1741 total_friction_force_.resize(nb_particles_tot,dimension);
1742 proportion_fa7_ok_UP2_.resize(nb_particles_tot,dimension);
1743 prop_P2_fluid_compo_.resize(nb_particles_tot);
1744 U_P2_moy_.resize(nb_particles_tot, dimension);
1745 total_surface_interf_.resize(nb_particles_tot);
1746 Nb_fa7_tot_par_compo_.resize(nb_particles_tot);
1747 total_heat_transfer_.resize(nb_particles_tot);
1748
1749 total_pressure_force_=0;
1750 total_friction_force_=0;
1751 proportion_fa7_ok_UP2_=0;
1752 prop_P2_fluid_compo_=0;
1753 U_P2_moy_=0;
1754 total_surface_interf_=0;
1755 Nb_fa7_tot_par_compo_=0;
1756}
1757
1758void Post_Processing_Hydrodynamic_Forces::resize_data_fa7(int nb_fa7)
1759{
1760 phase_indicator_function_P2_.resize(nb_fa7);
1761
1762 if (is_post_process_pressure_fa7_)
1763 {
1764 pressure_fa7_.resize(nb_fa7);
1765 pressure_fa7_=3e15;
1766 }
1767
1768 if (is_post_process_friction_force_fa7_)
1769 {
1770 friction_force_fa7_.resize(nb_fa7,dimension);
1771 friction_force_fa7_=3e15;
1772 }
1773
1774 if (is_post_process_pressure_force_fa7_)
1775 {
1776 pressure_force_fa7_.resize(nb_fa7,dimension);
1777 pressure_force_fa7_=3e15;
1778 }
1779
1780 if (is_post_process_stress_tensor_fa7_)
1781 {
1782 resize_sigma(nb_fa7);
1783 resize_gradU_P1(nb_fa7);
1784 resize_gradU_P2(nb_fa7);
1785 }
1786 U_P1_.resize(nb_fa7,dimension);
1787 U_P2_.resize(nb_fa7,dimension);
1788}
1789
1797
1798void Post_Processing_Hydrodynamic_Forces::compute_friction_force_complet_tensor(int nb_fa7,
1799 const Maillage_FT_Disc& mesh,
1800 const IntVect& compo_connexes_fa7,
1801 const ArrOfDouble& fa7_surface,
1802 const DoubleTab& tab_fa7_normal)
1803{
1804 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
1805 DoubleTab grad_U_P1(nb_fa7, dimension, dimension);
1806 DoubleTab grad_U_P2(nb_fa7, dimension, dimension);
1807 grad_U_P1=-1e15;
1808 grad_U_P2=-1e30;
1809 int interp_gradU_P1_ok=0;
1810 int interp_gradU_P2_ok=0;
1811 if (location_stress_tensor_== Location_stress_tensor::FACES_NORMALE_X)
1812 {
1813 interp_gradU_P1_ok=trilinear_interpolation_gradU_face(eq_ns.la_vitesse->valeurs(),
1814 coord_neighbor_fluid_fa7_gradU_1_, grad_U_P1);
1815 interp_gradU_P2_ok=trilinear_interpolation_gradU_face(eq_ns.la_vitesse->valeurs(),
1816 coord_neighbor_fluid_fa7_gradU_2_, grad_U_P2);
1817 }
1818 else if (location_stress_tensor_== Location_stress_tensor::ELEMENTS)
1819 {
1820 interp_gradU_P1_ok=trilinear_interpolation_gradU_elem(eq_ns.la_vitesse->valeurs(),
1821 coord_neighbor_fluid_fa7_gradU_1_,
1822 grad_U_P1);
1823 interp_gradU_P2_ok=trilinear_interpolation_gradU_elem(eq_ns.la_vitesse->valeurs(),
1824 coord_neighbor_fluid_fa7_gradU_2_,
1825 grad_U_P2);
1826 }
1827 if ( interp_gradU_P1_ok==1 && interp_gradU_P2_ok==1 )
1828 {
1829 DoubleTab grad_U_extrapole(nb_fa7, dimension, dimension);
1830 grad_U_extrapole=1e20;
1831
1832 for (int fa7=0; fa7<nb_fa7; fa7++)
1833 {
1834 fill_gradU_P1(fa7, grad_U_P1);
1835 fill_gradU_P2(fa7, grad_U_P2);
1836
1837 for (int i=0; i<dimension; i++)
1838 {
1839 for (int j=0; j<dimension; j++)
1840 {
1841 grad_U_extrapole(fa7,i,j) = grad_U_P2(fa7,i,j) -
1842 interpolation_distance_gradU_P2_*
1843 (grad_U_P2(fa7,i,j)-grad_U_P1(fa7,i,j))/
1844 (interpolation_distance_gradU_P2_-
1845 interpolation_distance_gradU_P1_);
1846 }
1847 }
1848 }
1849 for (int fa7=0; fa7<nb_fa7; fa7++)
1850 {
1851 int compo=compo_connexes_fa7(fa7);
1852 Matrice_Dense stress_tensor(dimension,dimension);
1853 for (int i=0; i<dimension; i++)
1854 {
1855 for (int j=0; j<dimension; j++)
1856 {
1857 stress_tensor(i,j) = (grad_U_extrapole(fa7,i,j) +
1858 grad_U_extrapole(fa7,j,i));
1859 }
1860 }
1861
1862 if (is_post_process_stress_tensor_fa7_)
1863 fill_sigma(fa7,stress_tensor);
1864
1865 DoubleTab la_normale_fa7_x_surface(dimension);
1866 for (int dim=0; dim<dimension; dim++) la_normale_fa7_x_surface(dim) =
1867 fa7_surface(fa7)*tab_fa7_normal(fa7,dim);
1868 DoubleVect friction_force_fa7=stress_tensor*la_normale_fa7_x_surface;
1869
1870 if (!mesh.facette_virtuelle(fa7))
1871 {
1872 if (is_post_process_friction_force_fa7_)
1873 {
1874 for (int dim=0; dim<dimension; dim++)
1875 friction_force_fa7_(fa7,dim)=friction_force_fa7(dim);
1876 }
1877 for (int dim=0; dim<dimension; dim++)
1878 total_friction_force_(compo,dim)+=friction_force_fa7(dim);
1879 }
1880 }
1881 }
1882}
1883
1884void Post_Processing_Hydrodynamic_Forces::compute_friction_force_projected_tensor(int nb_fa7,
1885 const Maillage_FT_Disc& mesh,
1886 const IntVect& compo_connexes_fa7,
1887 const ArrOfDouble& fa7_surface,
1888 const DoubleTab& tab_fa7_normal)
1889{
1890 Transport_Interfaces_FT_Disc& eq_transport = ptr_eq_transport_.valeur();
1891 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
1892 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, eq_ns.domaine_dis());
1893 const Fluide_Diphasique& two_phase_fluid = eq_ns.fluide_diphasique();
1894 const Domaine& domain = domain_vdf.domaine();
1895 const int id_fluid_phase=two_phase_fluid.get_id_fluid_phase();
1896 const double mu_f =two_phase_fluid.fluide_phase(id_fluid_phase).
1897 viscosite_dynamique().valeurs()(0, 0);
1898 const DoubleTab& gravity_center_fa7=mesh.get_gravity_center_fa7();
1899
1900 int interp_U_P1_ok=0;
1901 int interp_U_P2_ok=0;
1902 U_P1_.resize(nb_fa7, dimension);
1903 U_P2_.resize(nb_fa7, dimension);
1904
1905 DoubleTab U_P1_spherique(nb_fa7, dimension);
1906 DoubleTab U_P2_spherique(nb_fa7, dimension);
1907 DoubleTab U_cg_spherique(nb_fa7, dimension);
1908 DoubleTab Urr(nb_fa7);
1909 DoubleTab Uthetar(nb_fa7);
1910 DoubleTab Uphir(nb_fa7);
1911
1912 U_P1_=-1e15;
1913 U_P2_=-1e30;
1914 U_P1_spherique=-1e15;
1915 U_P2_spherique=-1e30;
1916 U_cg_spherique=-1e20;
1917 Urr=1e8;
1918 Uthetar=1e12;
1919 Uphir=1e15;
1920
1921 double theta=0;
1922 double phi=0;
1923 double distance_au_cg=0;
1924 const DoubleTab& positions_compo=eq_transport.get_particles_position();
1925 const DoubleTab& vitesses_compo = eq_transport.get_particles_velocity();
1926
1927 if (location_stress_tensor_== Location_stress_tensor::ELEMENTS)
1928 {
1929 // 1. Trilinear interpolation in cartesian coordinates
1930 interp_U_P1_ok=trilinear_interpolation_face(eq_ns.la_vitesse->valeurs(),
1931 coord_neighbor_fluid_fa7_gradU_1_,
1932 U_P1_);
1933 interp_U_P2_ok=trilinear_interpolation_face(eq_ns.la_vitesse->valeurs(),
1934 coord_neighbor_fluid_fa7_gradU_2_,
1935 U_P2_);
1936 }
1937 if ( interp_U_P1_ok && interp_U_P2_ok )
1938 {
1939 // 2. Change of reference frame to spherical coordinates
1940 for (int fa7=0; fa7<nb_fa7; fa7++)
1941 {
1942 int compo=compo_connexes_fa7(fa7);
1943 if (!mesh.facette_virtuelle(fa7))
1944 {
1945 Nb_fa7_tot_par_compo_(compo)++;
1946
1947 int cont=0; // if the velocity in P2 could not be computed
1948 // we do not compute the friction force for the fa7
1949 for (int dim=0; dim<dimension; dim++)
1950 {
1951 if (U_P2_(fa7,dim)>-1e10)
1952 {
1953 U_P2_moy_(compo,dim)+= U_P2_(fa7,dim); // compute the mean velocity in P2
1954 proportion_fa7_ok_UP2_(compo,dim)+=1;
1955 }
1956 else
1957 cont=1;
1958 }
1959 if (cont)
1960 continue;
1961
1962 DoubleVect distance_cg_vect(dimension);
1963 for (int i=0; i<dimension; i++) distance_cg_vect(i)=
1964 coord_neighbor_fluid_fa7_gradU_1_(fa7,i)-positions_compo(compo,i);
1965
1966 distance_au_cg=sqrt(local_carre_norme_vect(distance_cg_vect));
1967
1968 if (fabs((coord_neighbor_fluid_fa7_gradU_1_(fa7,2)-positions_compo(compo,2))/
1969 distance_au_cg)<=1)
1970 {
1971 theta=acos((coord_neighbor_fluid_fa7_gradU_1_(fa7,2)-
1972 positions_compo(compo,2))/distance_au_cg);
1973 }
1974 else if ((coord_neighbor_fluid_fa7_gradU_1_(fa7,2)-positions_compo(compo,2))/
1975 distance_au_cg>1)
1976 {
1977 theta=0;
1978 }
1979 else if ((coord_neighbor_fluid_fa7_gradU_1_(fa7,2)-positions_compo(compo,2))/
1980 distance_au_cg<-1)
1981 {
1982 theta=M_PI;
1983 }
1984
1985 if ((coord_neighbor_fluid_fa7_gradU_1_(fa7,0)-positions_compo(compo,0))>0 &&
1986 (coord_neighbor_fluid_fa7_gradU_1_(fa7,1)-positions_compo(compo,1))>=0)
1987 {
1988 phi=atan((coord_neighbor_fluid_fa7_gradU_1_(fa7,1)-positions_compo(compo,1))/
1989 (coord_neighbor_fluid_fa7_gradU_1_(fa7,0)-positions_compo(compo,0)));
1990 }
1991 else if ((coord_neighbor_fluid_fa7_gradU_1_(fa7,0)-positions_compo(compo,0))>0
1992 && (coord_neighbor_fluid_fa7_gradU_1_(fa7,1)-positions_compo(compo,1))<0)
1993 {
1994 phi=atan((coord_neighbor_fluid_fa7_gradU_1_(fa7,1)-positions_compo(compo,1))/
1995 (coord_neighbor_fluid_fa7_gradU_1_(fa7,0)-positions_compo(compo,0)))
1996 +2*M_PI;
1997 }
1998 else if ((coord_neighbor_fluid_fa7_gradU_1_(fa7,0)-positions_compo(compo,0))<0)
1999 {
2000 phi=atan((coord_neighbor_fluid_fa7_gradU_1_(fa7,1)-positions_compo(compo,1))/
2001 (coord_neighbor_fluid_fa7_gradU_1_(fa7,0)-positions_compo(compo,0)))
2002 +M_PI;
2003 }
2004 else if ((coord_neighbor_fluid_fa7_gradU_1_(fa7,0)-positions_compo(compo,0))==0
2005 && (coord_neighbor_fluid_fa7_gradU_1_(fa7,1)-positions_compo(compo,1))>0)
2006 {
2007 phi=M_PI/2.;
2008 }
2009 else if ((coord_neighbor_fluid_fa7_gradU_1_(fa7,0)-positions_compo(compo,0))==0
2010 && (coord_neighbor_fluid_fa7_gradU_1_(fa7,1)-positions_compo(compo,1))<0)
2011 {
2012 phi=3.*M_PI/2.;
2013 }
2014
2015 U_P1_spherique(fa7,0)=sin(theta)*cos(phi)*U_P1_(fa7,0)+sin(theta)*sin(phi)*
2016 U_P1_(fa7,1)+cos(theta)*U_P1_(fa7,2);
2017 U_P1_spherique(fa7,1)=cos(theta)*cos(phi)*U_P1_(fa7,0)+cos(theta)*sin(phi)*
2018 U_P1_(fa7,1)-sin(theta)*U_P1_(fa7,2);
2019 U_P1_spherique(fa7,2)=sin(phi)*U_P1_(fa7,0)+cos(phi)*U_P1_(fa7,1);
2020
2021 U_P2_spherique(fa7,0)=sin(theta)*cos(phi)*U_P2_(fa7,0)+sin(theta)*sin(phi)*
2022 U_P2_(fa7,1)+cos(theta)*U_P2_(fa7,2);
2023 U_P2_spherique(fa7,1)=cos(theta)*cos(phi)*U_P2_(fa7,0)+cos(theta)*sin(phi)*
2024 U_P2_(fa7,1)-sin(theta)*U_P2_(fa7,2);
2025 U_P2_spherique(fa7,2)=sin(phi)*U_P2_(fa7,0)+cos(phi)*U_P2_(fa7,1);
2026
2027 U_cg_spherique(fa7,0)=sin(theta)*cos(phi)*vitesses_compo(compo,0)+sin(theta)*
2028 sin(phi)*vitesses_compo(compo,1)+cos(theta)*vitesses_compo(compo,2);
2029 U_cg_spherique(fa7,1)=cos(theta)*cos(phi)*vitesses_compo(compo,0)+cos(theta)*
2030 sin(phi)*vitesses_compo(compo,1)-sin(theta)*vitesses_compo(compo,2);
2031 U_cg_spherique(fa7,2)=sin(phi)*vitesses_compo(compo,0)+cos(phi)*
2032 vitesses_compo(compo,1);
2033
2034 // we compute delta again because we need epsilon
2035 int elem_diph=domain.chercher_elements(gravity_center_fa7(fa7,0),
2036 gravity_center_fa7(fa7,1),
2037 gravity_center_fa7(fa7,2));
2038
2039 DoubleVect delta_i(dimension);
2040 int elem00=face_voisins_for_interp(elem_faces_for_interp(elem_diph, 0+dimension),1);
2041 int elem11=face_voisins_for_interp(elem_faces_for_interp(elem_diph, 1+dimension),1);
2042 int elem22=face_voisins_for_interp(elem_faces_for_interp(elem_diph, 2+dimension),1);
2043 int elem33=face_voisins_for_interp(elem_faces_for_interp(elem_diph, 2),0);
2044
2045 delta_i(0) = elem00>=0 ? fabs(domain_vdf.dist_elem(elem_diph, elem00, 0)) : -1e15;
2046 delta_i(1) = elem11>=0 ? fabs(domain_vdf.dist_elem(elem_diph, elem11, 1)) : -1e15;
2047
2048 if (tab_fa7_normal(fa7,2)>0)
2049 delta_i(2) = elem22>=0 ? fabs(domain_vdf.dist_elem(elem_diph, elem22, 2)) : -1e15;
2050 else
2051 delta_i(2) = elem33>=0 ? fabs(domain_vdf.dist_elem(elem_diph, elem33 , 2)) : -1e15;
2052
2053 double epsilon=0;
2054 // Interpolation distance varies with mesh refinement
2055 for (int dim=0; dim<dimension; dim++)
2056 epsilon+= fabs(delta_i(dim)*fabs(tab_fa7_normal(fa7,dim)));
2057
2058 // We compute the components of the friction force in spherical coordinates
2059 // afet simplification: ff=mu*(2*Urr, Uthetar, Uphir)
2060 Urr(fa7)=(-U_P2_spherique(fa7,0)+4.*U_P1_spherique(fa7,0)-3.*
2061 U_cg_spherique(fa7,0))/(2.*epsilon);
2062 Uthetar(fa7)=(-U_P2_spherique(fa7,1)+4.*U_P1_spherique(fa7,1)-3.*
2063 U_cg_spherique(fa7,1))/(2.*epsilon);
2064 Uphir(fa7)=(-U_P2_spherique(fa7,2)+4.*U_P1_spherique(fa7,2)-3.*
2065 U_cg_spherique(fa7,2))/(2.*epsilon);
2066
2067 DoubleVect ff(dimension);
2068 ff(0)=mu_f*fa7_surface(fa7)*(2.*sin(theta)*cos(phi)*Urr(fa7)+cos(theta)*
2069 cos(phi)*Uthetar(fa7)-sin(phi)*Uphir(fa7));
2070 ff(1)=mu_f*fa7_surface(fa7)*(2.*sin(theta)*sin(phi)*Urr(fa7)+cos(theta)*
2071 sin(phi)*Uthetar(fa7)+cos(phi)*Uphir(fa7));
2072 ff(2)=mu_f*fa7_surface(fa7)*(2.*cos(theta)*Urr(fa7)-sin(theta)*Uthetar(fa7));
2073
2074 for (int dim=0; dim<dimension; dim++)
2075 total_friction_force_(compo,dim)+=ff(dim);
2076
2077 if (is_post_process_friction_force_fa7_)
2078 {
2079 for (int dim=0; dim<dimension; dim++)
2080 friction_force_fa7_(fa7,dim)=ff(dim);
2081 }
2082 }
2083 }
2084 }
2085}
2086
2087void Post_Processing_Hydrodynamic_Forces::fill_absurd_value_coord_neighbor_fluid_fa7(int fa7)
2088{
2089 for (int dim=0; dim<dimension; dim++)
2090 {
2091 coord_neighbor_fluid_fa7_pressure_1_(fa7,dim)=-1e15;
2092 coord_neighbor_fluid_fa7_pressure_2_(fa7,dim)=-1e15;
2093 coord_neighbor_fluid_fa7_gradU_1_(fa7,dim)=-1e15;
2094 coord_neighbor_fluid_fa7_gradU_2_(fa7,dim)=-1e15;
2095 }
2096}
2097
2098void Post_Processing_Hydrodynamic_Forces::compute_pressure_force_trilinear_linear(int nb_fa7,
2099 const Maillage_FT_Disc& mesh,
2100 Convection_Diffusion_Temperature_FT_Disc::
2101 Thermal_correction_discretization_method
2102 thermal_correction_discretization_method,
2103 const IntVect& compo_connexes_fa7,
2104 const ArrOfDouble& fa7_surface,
2105 const DoubleTab& tab_fa7_normal)
2106{
2107 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
2108 DoubleTab pressure_P1(nb_fa7);
2109 DoubleTab pressure_P2(nb_fa7);
2110
2111 if (trilinear_interpolation_elem(eq_ns.la_pression->valeurs(),
2112 coord_neighbor_fluid_fa7_pressure_1_, pressure_P1, 0, thermal_correction_discretization_method) &&
2113 trilinear_interpolation_elem(eq_ns.la_pression->valeurs(),
2114 coord_neighbor_fluid_fa7_pressure_2_, pressure_P2,1,
2115 Convection_Diffusion_Temperature_FT_Disc::
2116 Thermal_correction_discretization_method::P1))
2117 {
2118 DoubleTab extrapolated_pressure(nb_fa7);
2119 prop_P2_fluid_compo_=0;
2120 for (int fa7=0; fa7<nb_fa7; fa7++)
2121 {
2122 if (!mesh.facette_virtuelle(fa7))
2123 {
2124 int compo = compo_connexes_fa7(fa7);
2125 if (phase_indicator_function_P2_(fa7)==1) prop_P2_fluid_compo_(compo)+=1;
2126
2127 if (pressure_P2(fa7)<-1e10)
2128 {
2129 if (is_post_process_pressure_fa7_) pressure_fa7_(fa7)=1e15;
2130 extrapolated_pressure(fa7)=1e15;
2131 continue;
2132 }
2133
2134 // 3*pressure_P2(compo,fa7)-2*pressure_P1(compo,fa7); //
2135 // If one cannot interpolate in P1 and P2, then one do not compute the pressure force
2136 extrapolated_pressure(fa7) = pressure_P2(fa7)-interpolation_distance_pressure_P2_*
2137 (pressure_P2(fa7)-pressure_P1(fa7))/(interpolation_distance_pressure_P2_-
2138 interpolation_distance_pressure_P1_);
2139 if (is_post_process_pressure_fa7_)
2140 pressure_fa7_(fa7)=extrapolated_pressure(fa7);
2141 }
2142 else
2143 {
2144 if (is_post_process_pressure_fa7_) pressure_fa7_(fa7)=1e15;
2145 extrapolated_pressure(fa7)=-1e15;
2146 }
2147 }
2148 for (int fa7=0; fa7<nb_fa7; fa7++)
2149 {
2150 int compo = compo_connexes_fa7(fa7);
2151 if (extrapolated_pressure(fa7)>1e10)
2152 continue;
2153 double coeff=-extrapolated_pressure(fa7)*fa7_surface(fa7);
2154 DoubleVect pressure_force_fa7(dimension);
2155 for (int dim=0; dim<dimension; dim++)
2156 pressure_force_fa7(dim)=coeff*tab_fa7_normal(fa7,dim);
2157 if (!mesh.facette_virtuelle(fa7))
2158 {
2159 total_surface_interf_(compo)+=fa7_surface(fa7);
2160 if (is_post_process_pressure_force_fa7_)
2161 {
2162 for (int dim=0; dim<dimension; dim++)
2163 pressure_force_fa7_(fa7,dim)=pressure_force_fa7(dim);
2164 }
2165 for (int dim=0; dim<dimension; dim++)
2166 total_pressure_force_(compo,dim)+=pressure_force_fa7(dim);
2167 }
2168 }
2169 }
2170}
2171
2172void Post_Processing_Hydrodynamic_Forces::compute_neighbors_coordinates_fluid_fa7(const int nb_fa7,
2173 const int is_discr_elem_diph,
2174 const DoubleTab& gravity_center_fa7,
2175 const Maillage_FT_Disc& mesh,
2176 const DoubleTab& tab_fa7_normal,
2177 const IntTab& particles_eulerian_id_number)
2178{
2179 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
2180 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, eq_ns.domaine_dis());
2181 const Domaine& domain = domain_vdf.domaine();
2182
2183 for (int fa7 =0 ; fa7<nb_fa7 ; fa7++)
2184 {
2185 if (!mesh.facette_virtuelle(fa7))
2186 {
2187 DoubleVect fa7_normal(dimension);
2188 int elem_diph=domain.chercher_elements(gravity_center_fa7(fa7,0),
2189 gravity_center_fa7(fa7,1),
2190 gravity_center_fa7(fa7,2));
2191 if (is_discr_elem_diph)
2192 {
2193 list_elem_diph_(fa7,0) = elem_diph;
2194 list_elem_diph_(fa7,1) = particles_eulerian_id_number(elem_diph);
2195 }
2196 if (elem_diph>=0)
2197 {
2198 DoubleVect delta_i(dimension);
2199 // One compute eulerian mesh thicknesses in which are located the lagrangian facets
2200 // If we have access, one compute the thickness outside of the particle
2201 // Otherwise, one compute the thickness inside of the particle
2202 // It just comes to the selection of the mesh juxtaposed to the two-phase cell.
2203 int acces=1;
2204 for (int dim=0; dim<dimension; dim++)
2205 {
2206 int elem_haut=face_voisins_for_interp(
2207 elem_faces_for_interp(elem_diph, dim+dimension),1);
2208 int elem_bas=face_voisins_for_interp(
2209 elem_faces_for_interp(elem_diph, dim),0);
2210 if (elem_bas<0 && elem_haut<0)
2211 acces=0;
2212 if (tab_fa7_normal(fa7,dim)>0)
2213 {
2214 delta_i(dim) = (elem_haut>=0) ? fabs(domain_vdf.dist_elem(elem_diph,elem_haut, dim)) :
2215 fabs(domain_vdf.dist_elem(elem_diph,elem_bas, dim));
2216 }
2217 else
2218 {
2219 delta_i(dim) = (elem_bas>=0) ? fabs(domain_vdf.dist_elem(elem_diph,elem_bas, dim)) :
2220 fabs(domain_vdf.dist_elem(elem_diph,elem_haut, dim)) ;
2221 }
2222 }
2223 if (acces==0)
2224 fill_absurd_value_coord_neighbor_fluid_fa7(fa7);
2225 else
2226 {
2227 double epsilon=0;
2228 for (int dim=0; dim<dimension; dim++)
2229 {
2230 // Interpolation distance varies with mesh refinement
2231 epsilon+= fabs(delta_i(dim)*fabs(tab_fa7_normal(fa7,dim)));
2232 }
2233 for (int dim=0; dim<dimension; dim++)
2234 {
2235 fa7_normal(dim)=tab_fa7_normal(fa7,dim);
2236 coord_neighbor_fluid_fa7_pressure_1_(fa7,dim)=gravity_center_fa7(fa7,dim)+
2237 interpolation_distance_pressure_P1_*epsilon*fa7_normal(dim);
2238 coord_neighbor_fluid_fa7_pressure_2_(fa7,dim)=gravity_center_fa7(fa7,dim)+
2239 interpolation_distance_pressure_P2_*epsilon*fa7_normal(dim);
2240 coord_neighbor_fluid_fa7_gradU_1_(fa7,dim)=gravity_center_fa7(fa7,dim)+
2241 interpolation_distance_pressure_P1_*epsilon*fa7_normal(dim);
2242 coord_neighbor_fluid_fa7_gradU_2_(fa7,dim)=gravity_center_fa7(fa7,dim)+
2243 interpolation_distance_pressure_P2_*epsilon*fa7_normal(dim);
2244 }
2245 }
2246 }
2247 else
2248 fill_absurd_value_coord_neighbor_fluid_fa7(fa7);
2249 }
2250 }
2251}
2252
2253void Post_Processing_Hydrodynamic_Forces::compute_U_P2_moy(const int nb_particles_tot)
2254{
2255 for (int compo=0; compo<nb_particles_tot; compo++)
2256 {
2257 for (int dim=0; dim<dimension; dim++)
2258 {
2259 if (proportion_fa7_ok_UP2_(compo,dim)>0)
2260 U_P2_moy_(compo,dim)/=proportion_fa7_ok_UP2_(compo,dim);
2261 }
2262 }
2263}
2264
2265void Post_Processing_Hydrodynamic_Forces::compute_proportion_fa7_ok_and_is_fluid_P2
2266(const int nb_particles_tot)
2267{
2268 for (int compo=0; compo<nb_particles_tot; compo++)
2269 {
2270 for (int dim=0; dim<dimension; dim++)
2271 {
2272 if (Nb_fa7_tot_par_compo_(compo)>0)
2273 proportion_fa7_ok_UP2_(compo,dim)/=Nb_fa7_tot_par_compo_(compo);
2274 }
2275 if (Nb_fa7_tot_par_compo_(compo)>0)
2276 prop_P2_fluid_compo_(compo)/=Nb_fa7_tot_par_compo_(compo);
2277
2278 }
2279}
2280
2281void Post_Processing_Hydrodynamic_Forces::compute_heat_transfer()
2282{
2283 if(!flag_heat_transfer_computation_)
2284 {
2285 Cerr << "Post_Processing_Hydrodynamic_Forces::compute_heat_transfer" << finl;
2286 const Navier_Stokes_FT_Disc& eq_ns = ptr_eq_ns_.valeur();
2287 Transport_Interfaces_FT_Disc& eq_transport = ptr_eq_transport_.valeur();
2288 Convection_Diffusion_Temperature_FT_Disc& eq_temp = ptr_eq_temp_.valeur();
2289 const DoubleTab& temperature = eq_temp.inconnue().valeurs();
2290 const Maillage_FT_Disc& mesh = eq_transport.maillage_interface();
2291 const Domaine_VDF& domain_vdf = ref_cast(Domaine_VDF, eq_ns.domaine_dis());
2292 const Domaine& domain = domain_vdf.domaine();
2293 const Fluide_Diphasique& mon_fluide = eq_ns.fluide_diphasique();
2294 double lambda_f=mon_fluide.fluide_phase(1).conductivite().valeurs()(0, 0);
2295 const int nb_fa7 = mesh.nb_facettes();
2296
2297 IntVect compo_connexes_fa7(nb_fa7);
2298 int n = search_connex_components_local_FT(mesh, compo_connexes_fa7);
2299 int nb_compo_tot=compute_global_connex_components_FT(mesh, compo_connexes_fa7, n);
2300
2301 if (total_heat_transfer_.size_array() != nb_compo_tot)
2302 {
2303 total_heat_transfer_.resize(nb_compo_tot);
2304 total_heat_transfer_=0;
2305 }
2306
2307 if (nb_fa7>0)
2308 {
2309 const ArrOfDouble& les_surfaces_fa7 = mesh.get_update_surface_facettes();
2310 const DoubleTab& les_normales_fa7 = mesh.get_update_normale_facettes();
2311 heat_transfer_fa7_.resize(nb_fa7);
2312 heat_transfer_fa7_=1e15;
2313
2314 const DoubleTab& les_cg_fa7=mesh.get_gravity_center_fa7();
2315 coord_neighbor_fluid_fa7_temp_1_.resize(nb_fa7,dimension);
2316 coord_neighbor_fluid_fa7_temp_2_.resize(nb_fa7,dimension);
2317
2318 for (int fa7 =0 ; fa7<nb_fa7 ; fa7++)
2319 {
2320 if (!mesh.facette_virtuelle(fa7))
2321 {
2322 DoubleVect normale_fa7(dimension);
2323 int elem_diph=domain.chercher_elements(les_cg_fa7(fa7,0),
2324 les_cg_fa7(fa7,1),les_cg_fa7(fa7,2));
2325 DoubleVect delta_i(dimension);
2326 for (int dim=0; dim<dimension; dim++)
2327 {
2328 int elem_haut=face_voisins_for_interp(elem_faces_for_interp(elem_diph,
2329 dim+dimension),1);
2330 int elem_bas=face_voisins_for_interp(elem_faces_for_interp(elem_diph, dim),0);
2331 if (les_normales_fa7(fa7,dim)>0)
2332 delta_i(dim) = (elem_haut>=0) ? fabs(domain_vdf.dist_elem(elem_diph,
2333 elem_haut, dim)) : fabs(domain_vdf.dist_elem(elem_diph,elem_bas, dim));
2334 else delta_i(dim) = (elem_bas>=0) ? fabs(domain_vdf.dist_elem(elem_diph,
2335 elem_bas, dim)) : fabs(domain_vdf.dist_elem(elem_diph,elem_haut, dim));
2336 }
2337 double epsilon=0;
2338 for (int dim=0; dim<dimension; dim++)
2339 epsilon+= fabs(delta_i(dim)*fabs(les_normales_fa7(fa7,dim)));
2340 for (int dim=0; dim<dimension; dim++)
2341 {
2342 normale_fa7(dim)=les_normales_fa7(fa7,dim);
2343 coord_neighbor_fluid_fa7_temp_1_(fa7,dim)=les_cg_fa7(fa7,dim)+
2344 interpolation_distance_temperature_P1_*epsilon*normale_fa7(dim);
2345 coord_neighbor_fluid_fa7_temp_2_(fa7,dim)=les_cg_fa7(fa7,dim)+
2346 interpolation_distance_temperature_P2_*epsilon*normale_fa7(dim);
2347 }
2348 }
2349 }
2350
2351 DoubleTab temp_P1(nb_fa7);
2352 DoubleTab temp_P2(nb_fa7);
2353
2354 int interp_T_P1_ok=trilinear_interpolation_elem(temperature,
2355 coord_neighbor_fluid_fa7_temp_1_, temp_P1);
2356 int interp_T_P2_ok=trilinear_interpolation_elem(temperature,
2357 coord_neighbor_fluid_fa7_temp_2_, temp_P2);
2358 if (interp_T_P1_ok && interp_T_P2_ok)
2359 {
2360 for (int fa7=0; fa7<nb_fa7; fa7++)
2361 {
2362 int compo=compo_connexes_fa7(fa7);
2363 if (!mesh.facette_virtuelle(fa7))
2364 {
2365 int elem_diph=domain.chercher_elements(les_cg_fa7(fa7,0),
2366 les_cg_fa7(fa7,1),les_cg_fa7(fa7,2));
2367 DoubleVect delta_i(dimension);
2368 int elem00=face_voisins_for_interp(elem_faces_for_interp(elem_diph, 0+dimension),1);
2369 int elem11=face_voisins_for_interp(elem_faces_for_interp(elem_diph, 1+dimension),1);
2370 int elem22=face_voisins_for_interp(elem_faces_for_interp(elem_diph, 2+dimension),1);
2371 int elem33=face_voisins_for_interp(elem_faces_for_interp(elem_diph, 2),0);
2372
2373 delta_i(0) = elem00>=0 ? fabs(domain_vdf.dist_elem(elem_diph, elem00, 0)) : -1e15;
2374 delta_i(1) = elem11>=0 ? fabs(domain_vdf.dist_elem(elem_diph, elem11, 1)) : -1e15;
2375
2376 if (les_normales_fa7(fa7,2)>0)
2377 delta_i(2) = elem22>=0 ? fabs(domain_vdf.dist_elem(elem_diph, elem22, 2)) : -1e15;
2378 else
2379 delta_i(2) = elem33>=0 ? fabs(domain_vdf.dist_elem(elem_diph, elem33 , 2)) : -1e15;
2380
2381 double epsilon=0;
2382 for (int dim=0; dim<dimension; dim++)
2383 epsilon+= fabs(delta_i(dim)*fabs(les_normales_fa7(fa7,dim))); // la distance d'interpolation varie en fonction du raffinement du maillage
2384 heat_transfer_fa7_(fa7)=lambda_f*(-temp_P2(fa7)+4.*temp_P1(fa7)-3.*eq_temp.get_tsat_constant())/(2.*epsilon)*les_surfaces_fa7(fa7); // schema decentre avant d'ordre 2
2385 total_heat_transfer_(compo)+=heat_transfer_fa7_(fa7);
2386 }
2387 }
2388 }
2389 }
2390 mp_sum_for_each_item(total_heat_transfer_);
2391 raise_the_flag_heat_transfer();
2392 }
2393}
2394
Champ_Don_base::valeurs
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
Definition Champ_Don_base.h:49
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
Definition Champ_Inc_base.h:77
Champ_Proto::valeurs
virtual DoubleTab & valeurs()=0
Convection_Diffusion_Temperature_FT_Disc
Definition Convection_Diffusion_Temperature_FT_Disc.h:30
Convection_Diffusion_Temperature_FT_Disc::get_tsat_constant
const double & get_tsat_constant() const
Definition Convection_Diffusion_Temperature_FT_Disc.cpp:1943
Convection_Diffusion_Temperature_FT_Disc::Thermal_correction_discretization_method
Thermal_correction_discretization_method
Definition Convection_Diffusion_Temperature_FT_Disc.h:98
Convection_Diffusion_Temperature_FT_Disc::Thermal_correction_discretization_method::P1
@ P1
Definition Convection_Diffusion_Temperature_FT_Disc.h:98
Convection_Diffusion_Temperature::inconnue
const Champ_Inc_base & inconnue() const override
Definition Convection_Diffusion_Temperature.h:34
Domaine_32_64::chercher_elements
SmallArrOfTID_t & chercher_elements(const DoubleTab &pos, SmallArrOfTID_t &elem, int reel=0) const
Recherche des elements contenant les points dont les coordonnees sont specifiees.
Definition Domaine.cpp:405
Domaine_VDF
class Domaine_VDF
Definition Domaine_VDF.h:64
Domaine_VDF::dist_face
double dist_face(int, int, int k) const
Definition Domaine_VDF.h:308
Domaine_VDF::dist_elem
double dist_elem(int, int, int) const
Definition Domaine_VDF.h:610
Domaine_VF
class Domaine_VF
Definition Domaine_VF.h:44
Domaine_VF::xv
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
Domaine_VF::elem_faces
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
Definition Domaine_VF.h:543
Domaine_VF::xp
double xp(int num_elem, int k) const
Definition Domaine_VF.h:77
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Equation_base::domaine_dis
Domaine_dis_base & domaine_dis()
Renvoie le domaine discretise associe a l'equation.
Definition Equation_base.cpp:92
Fluide_Diphasique
Definition Fluide_Diphasique.h:26
Fluide_Diphasique::fluide_phase
const Fluide_Incompressible & fluide_phase(int la_phase) const
Definition Fluide_Diphasique.cpp:100
Fluide_Diphasique::get_id_fluid_phase
const int & get_id_fluid_phase() const
Definition Fluide_Diphasique.h:55
Maillage_FT_Disc
: class Maillage_FT_Disc Cette classe decrit un maillage:
Definition Maillage_FT_Disc.h:47
Maillage_FT_Disc::sommets
const DoubleTab & sommets() const
renvoie le tableau des sommets (reels et virtuels) dimension(0) = nombre de sommets,
Definition Maillage_FT_Disc.cpp:1816
Maillage_FT_Disc::nb_facettes
int nb_facettes() const
renvoie le nombre de facettes (reelles et virtuelles) (egal a facettes().
Definition Maillage_FT_Disc.cpp:1852
Maillage_FT_Disc::facette_virtuelle
int facette_virtuelle(int i) const
Renvoie 0 si la facette m'appartient, 1 sinon.
Definition Maillage_FT_Disc.h:528
Maillage_FT_Disc::get_update_normale_facettes
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...
Definition Maillage_FT_Disc.cpp:1586
Maillage_FT_Disc::get_gravity_center_fa7
const DoubleTab & get_gravity_center_fa7() const
Definition Maillage_FT_Disc.h:234
Maillage_FT_Disc::get_update_surface_facettes
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...
Definition Maillage_FT_Disc.cpp:1563
Maillage_FT_Disc::sommet_elem
const ArrOfInt & sommet_elem() const
pour postraitement, renvoie sommet_elem_
Definition Maillage_FT_Disc.cpp:2137
Matrice_Dense
Definition Matrice_Dense.h:25
Milieu_base::conductivite
virtual const Champ_Don_base & conductivite() const
Renvoie la conductivite du milieu.
Definition Milieu_base.cpp:847
Motcle
Une chaine de caractere (Nom) en majuscules.
Definition Motcle.h:26
Motcles
Un tableau d'objets de la classe Motcle.
Definition Motcle.h:63
Motcles::search
int search(const Motcle &t) const
Definition Motcle.cpp:321
Navier_Stokes_FT_Disc
Definition Navier_Stokes_FT_Disc.h:31
Navier_Stokes_FT_Disc::fluide_diphasique
virtual const Fluide_Diphasique & fluide_diphasique() const
Definition Navier_Stokes_FT_Disc.cpp:647
Navier_Stokes_FT_Disc::get_particles_eulerian_id_number
const IntTab & get_particles_eulerian_id_number() const
Definition Navier_Stokes_FT_Disc.h:79
Navier_Stokes_std::la_pression
la_pression
Definition Navier_Stokes_std.h:174
Objet_U_With_Params
Inherits from Objet_U, adds the very common method set_param for the Objet_U hierarchy.
Definition Objet_U_With_Params.h:44
Objet_U::set_param
virtual void set_param(Param &) const
Definition Objet_U.h:135
Objet_U::Entree
friend class Entree
Definition Objet_U.h:76
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Objet_U::que_suis_je
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Param
Helper class to factorize the readOn method of Objet_U classes.
Definition Param.h:112
Param::ajouter_flag
void ajouter_flag(const char *keyword, const bool *value)
Register a boolean flag whose mere presence switches it to true.
Definition Param.cpp:474
Param::ajouter
void ajouter(const char *keyword, const int *value, Param::Nature nat=Param::OPTIONAL)
Register an integer parameter.
Definition Param.cpp:364
Param::REQUIRED
@ REQUIRED
Definition Param.h:115
Param::ajouter_non_std
void ajouter_non_std(const char *keyword, const Objet_U *value, Param::Nature nat=Param::OPTIONAL)
Register a keyword handled by Objet_U::lire_motcle_non_standard.
Definition Param.cpp:489
Post_Processing_Hydrodynamic_Forces
Definition Post_Processing_Hydrodynamic_Forces.h:30
Post_Processing_Hydrodynamic_Forces::list_elem_P1_all_
IntTab list_elem_P1_all_
Definition Post_Processing_Hydrodynamic_Forces.h:237
Post_Processing_Hydrodynamic_Forces::coord_neighbor_fluid_fa7_gradU_2_
DoubleTab coord_neighbor_fluid_fa7_gradU_2_
Definition Post_Processing_Hydrodynamic_Forces.h:286
Post_Processing_Hydrodynamic_Forces::sigma_zz_fa7_
DoubleTab sigma_zz_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:261
Post_Processing_Hydrodynamic_Forces::lire_motcle_non_standard
int lire_motcle_non_standard(const Motcle &, Entree &) override
Lecture des parametres de type non simple d'un objet_U a partir d'un flot d'entree.
Definition Post_Processing_Hydrodynamic_Forces.cpp:63
Post_Processing_Hydrodynamic_Forces::raise_the_flag
void raise_the_flag()
Definition Post_Processing_Hydrodynamic_Forces.h:45
Post_Processing_Hydrodynamic_Forces::raise_the_flag_heat_transfer
void raise_the_flag_heat_transfer()
Definition Post_Processing_Hydrodynamic_Forces.h:47
Post_Processing_Hydrodynamic_Forces::dWdx_P1_
DoubleTab dWdx_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:269
Post_Processing_Hydrodynamic_Forces::total_pressure_force_
DoubleTab total_pressure_force_
Definition Post_Processing_Hydrodynamic_Forces.h:243
Post_Processing_Hydrodynamic_Forces::sigma_yy_fa7_
DoubleTab sigma_yy_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:257
Post_Processing_Hydrodynamic_Forces::fill_absurd_value_coord_neighbor_fluid_fa7
void fill_absurd_value_coord_neighbor_fluid_fa7(int fa7)
Definition Post_Processing_Hydrodynamic_Forces.cpp:2087
Post_Processing_Hydrodynamic_Forces::resize_sigma
virtual void resize_sigma(int nb_fa7)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1660
Post_Processing_Hydrodynamic_Forces::sigma_zy_fa7_
DoubleTab sigma_zy_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:260
Post_Processing_Hydrodynamic_Forces::list_elem_P1_
IntTab list_elem_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:238
Post_Processing_Hydrodynamic_Forces::find_neighboring_faces
void find_neighboring_faces(DoubleVect &coord_elem_interp, IntVect &neighboring_faces, int orientation) const
Definition Post_Processing_Hydrodynamic_Forces.cpp:1386
Post_Processing_Hydrodynamic_Forces::interpolation_distance_pressure_P2_
double interpolation_distance_pressure_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:227
Post_Processing_Hydrodynamic_Forces::U_P1_
DoubleTab U_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:250
Post_Processing_Hydrodynamic_Forces::sigma_xx_fa7_
DoubleTab sigma_xx_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:253
Post_Processing_Hydrodynamic_Forces::compute_friction_force_projected_tensor
virtual void compute_friction_force_projected_tensor(int nb_fa7, const Maillage_FT_Disc &mesh, const IntVect &compo_connexes_fa7, const ArrOfDouble &fa7_surface, const DoubleTab &tab_fa7_normal)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1884
Post_Processing_Hydrodynamic_Forces::compute_U_P2_moy
void compute_U_P2_moy(const int nb_particles_tot)
Definition Post_Processing_Hydrodynamic_Forces.cpp:2253
Post_Processing_Hydrodynamic_Forces::is_compute_heat_transfer_
bool is_compute_heat_transfer_
Definition Post_Processing_Hydrodynamic_Forces.h:218
Post_Processing_Hydrodynamic_Forces::is_post_process_pressure_force_fa7_
bool is_post_process_pressure_force_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:214
Post_Processing_Hydrodynamic_Forces::friction_force_fa7_
DoubleTab friction_force_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:247
Post_Processing_Hydrodynamic_Forces::compute_heat_transfer
void compute_heat_transfer()
Definition Post_Processing_Hydrodynamic_Forces.cpp:2281
Post_Processing_Hydrodynamic_Forces::heat_transfer_fa7_
DoubleTab heat_transfer_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:248
Post_Processing_Hydrodynamic_Forces::sigma_yx_fa7_
DoubleTab sigma_yx_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:256
Post_Processing_Hydrodynamic_Forces::U_P2_
DoubleTab U_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:251
Post_Processing_Hydrodynamic_Forces::dUdz_P1_
DoubleTab dUdz_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:265
Post_Processing_Hydrodynamic_Forces::coord_neighbor_fluid_fa7_pressure_1_
DoubleTab coord_neighbor_fluid_fa7_pressure_1_
Definition Post_Processing_Hydrodynamic_Forces.h:283
Post_Processing_Hydrodynamic_Forces::dUdy_P1_
DoubleTab dUdy_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:264
Post_Processing_Hydrodynamic_Forces::interpolation_distance_pressure_P1_
double interpolation_distance_pressure_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:226
Post_Processing_Hydrodynamic_Forces::coord_neighbor_fluid_fa7_temp_2_
DoubleTab coord_neighbor_fluid_fa7_temp_2_
Definition Post_Processing_Hydrodynamic_Forces.h:288
Post_Processing_Hydrodynamic_Forces::sigma_xy_fa7_
DoubleTab sigma_xy_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:254
Post_Processing_Hydrodynamic_Forces::compute_neighbors_coordinates_fluid_fa7
virtual void compute_neighbors_coordinates_fluid_fa7(const int nb_fa7, const int is_discr_elem_diph, const DoubleTab &gravity_center_fa7, const Maillage_FT_Disc &mesh, const DoubleTab &tab_fa7_normal, const IntTab &particles_eulerian_id_number)
Definition Post_Processing_Hydrodynamic_Forces.cpp:2172
Post_Processing_Hydrodynamic_Forces::compute_proportion_fa7_ok_and_is_fluid_P2
void compute_proportion_fa7_ok_and_is_fluid_P2(const int nb_particles_tot)
Definition Post_Processing_Hydrodynamic_Forces.cpp:2266
Post_Processing_Hydrodynamic_Forces::dWdx_P2_
DoubleTab dWdx_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:279
Post_Processing_Hydrodynamic_Forces::flag_heat_transfer_computation_
int flag_heat_transfer_computation_
Definition Post_Processing_Hydrodynamic_Forces.h:222
Post_Processing_Hydrodynamic_Forces::is_post_process_friction_force_fa7_
bool is_post_process_friction_force_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:215
Post_Processing_Hydrodynamic_Forces::dWdz_P1_
DoubleTab dWdz_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:271
Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_elem
int trilinear_interpolation_elem(const DoubleTab &valeurs_champ, DoubleTab &coord, DoubleTab &resu)
Definition Post_Processing_Hydrodynamic_Forces.cpp:255
Post_Processing_Hydrodynamic_Forces::dUdy_P2_
DoubleTab dUdy_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:274
Post_Processing_Hydrodynamic_Forces::fill_gradU_P2
void fill_gradU_P2(int fa7, DoubleTab gradU_P2)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1712
Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_face_sommets
int trilinear_interpolation_face_sommets(const DoubleTab &valeurs_champ, DoubleTab &coord, DoubleTab &resu)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1189
Post_Processing_Hydrodynamic_Forces::is_compute_forces_
bool is_compute_forces_
Definition Post_Processing_Hydrodynamic_Forces.h:219
Post_Processing_Hydrodynamic_Forces::dVdy_P2_
DoubleTab dVdy_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:277
Post_Processing_Hydrodynamic_Forces::dVdx_P2_
DoubleTab dVdx_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:276
Post_Processing_Hydrodynamic_Forces::dVdy_P1_
DoubleTab dVdy_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:267
Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_gradU_elem
int trilinear_interpolation_gradU_elem(const DoubleTab &valeurs_champ, DoubleTab &coord, DoubleTab &resu)
Definition Post_Processing_Hydrodynamic_Forces.cpp:985
Post_Processing_Hydrodynamic_Forces::sigma_yz_fa7_
DoubleTab sigma_yz_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:258
Post_Processing_Hydrodynamic_Forces::resize_gradU_P2
virtual void resize_gradU_P2(int nb_fa7)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1686
Post_Processing_Hydrodynamic_Forces::elem_faces_for_interp
int elem_faces_for_interp(int num_face, int i) const
Definition Post_Processing_Hydrodynamic_Forces.cpp:230
Post_Processing_Hydrodynamic_Forces::dWdy_P2_
DoubleTab dWdy_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:280
Post_Processing_Hydrodynamic_Forces::resize_data_fa7
virtual void resize_data_fa7(int nb_fa7)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1758
Post_Processing_Hydrodynamic_Forces::Method_pressure_force_computation::TRILINEAR_LINEAR
@ TRILINEAR_LINEAR
Definition Post_Processing_Hydrodynamic_Forces.h:112
Post_Processing_Hydrodynamic_Forces::dUdz_P2_
DoubleTab dUdz_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:275
Post_Processing_Hydrodynamic_Forces::fill_sigma
void fill_sigma(int fa7, Matrice_Dense stress_tensor)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1725
Post_Processing_Hydrodynamic_Forces::find_neighboring_faces_xyz
void find_neighboring_faces_xyz(DoubleVect &coord_elem_interp, IntTab &neighboring_faces) const
Definition Post_Processing_Hydrodynamic_Forces.cpp:1572
Post_Processing_Hydrodynamic_Forces::total_surface_interf_
DoubleVect total_surface_interf_
Definition Post_Processing_Hydrodynamic_Forces.h:233
Post_Processing_Hydrodynamic_Forces::sigma_xz_fa7_
DoubleTab sigma_xz_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:255
Post_Processing_Hydrodynamic_Forces::coord_neighbor_fluid_fa7_pressure_2_
DoubleTab coord_neighbor_fluid_fa7_pressure_2_
Definition Post_Processing_Hydrodynamic_Forces.h:284
Post_Processing_Hydrodynamic_Forces::dWdy_P1_
DoubleTab dWdy_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:270
Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_face
int trilinear_interpolation_face(const DoubleTab &valeurs_champ, DoubleTab &coord, DoubleTab &resu)
Definition Post_Processing_Hydrodynamic_Forces.cpp:406
Post_Processing_Hydrodynamic_Forces::Nb_fa7_tot_par_compo_
IntTab Nb_fa7_tot_par_compo_
Definition Post_Processing_Hydrodynamic_Forces.h:239
Post_Processing_Hydrodynamic_Forces::compute_hydrodynamic_forces
virtual void compute_hydrodynamic_forces()
Definition Post_Processing_Hydrodynamic_Forces.cpp:144
Post_Processing_Hydrodynamic_Forces::compute_pressure_force_trilinear_linear
virtual void compute_pressure_force_trilinear_linear(int nb_fa7, const Maillage_FT_Disc &mesh, Convection_Diffusion_Temperature_FT_Disc::Thermal_correction_discretization_method thermal_correction_discretization_method, const IntVect &compo_connexes_fa7, const ArrOfDouble &fa7_surface, const DoubleTab &tab_fa7_normal)
Definition Post_Processing_Hydrodynamic_Forces.cpp:2098
Post_Processing_Hydrodynamic_Forces::Location_stress_tensor::ELEMENTS
@ ELEMENTS
Definition Post_Processing_Hydrodynamic_Forces.h:113
Post_Processing_Hydrodynamic_Forces::Location_stress_tensor::FACES_NORMALE_X
@ FACES_NORMALE_X
Definition Post_Processing_Hydrodynamic_Forces.h:113
Post_Processing_Hydrodynamic_Forces::is_post_process_stress_tensor_fa7_
bool is_post_process_stress_tensor_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:216
Post_Processing_Hydrodynamic_Forces::dUdx_P2_
DoubleTab dUdx_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:273
Post_Processing_Hydrodynamic_Forces::location_stress_tensor_
Location_stress_tensor location_stress_tensor_
Definition Post_Processing_Hydrodynamic_Forces.h:179
Post_Processing_Hydrodynamic_Forces::total_friction_force_
DoubleTab total_friction_force_
Definition Post_Processing_Hydrodynamic_Forces.h:244
Post_Processing_Hydrodynamic_Forces::resize_gradU_P1
virtual void resize_gradU_P1(int nb_fa7)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1673
Post_Processing_Hydrodynamic_Forces::OBS_PTR
OBS_PTR(Transport_Interfaces_FT_Disc) ptr_eq_transport_
Post_Processing_Hydrodynamic_Forces::pressure_fa7_
DoubleVect pressure_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:235
Post_Processing_Hydrodynamic_Forces::is_compute_stokes_theoretical_forces_
bool is_compute_stokes_theoretical_forces_
Definition Post_Processing_Hydrodynamic_Forces.h:213
Post_Processing_Hydrodynamic_Forces::interpolation_distance_temperature_P2_
double interpolation_distance_temperature_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:225
Post_Processing_Hydrodynamic_Forces::dUdx_P1_
DoubleTab dUdx_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:263
Post_Processing_Hydrodynamic_Forces::U_P2_moy_
DoubleTab U_P2_moy_
Definition Post_Processing_Hydrodynamic_Forces.h:249
Post_Processing_Hydrodynamic_Forces::flag_force_computation_
int flag_force_computation_
Definition Post_Processing_Hydrodynamic_Forces.h:221
Post_Processing_Hydrodynamic_Forces::compute_friction_force_complet_tensor
virtual void compute_friction_force_complet_tensor(int nb_fa7, const Maillage_FT_Disc &mesh, const IntVect &compo_connexes_fa7, const ArrOfDouble &fa7_surface, const DoubleTab &tab_fa7_normal)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1798
Post_Processing_Hydrodynamic_Forces::phase_indicator_function_P2_
DoubleVect phase_indicator_function_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:232
Post_Processing_Hydrodynamic_Forces::resize_and_init_tables
virtual void resize_and_init_tables(int nb_compo_tot)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1738
Post_Processing_Hydrodynamic_Forces::find_neighboring_elements
double find_neighboring_elements(DoubleVect &coord_elem_interp, IntVect &neighboring_elements, const int sauv_list_P1=0, const int num_fa7=-1)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1271
Post_Processing_Hydrodynamic_Forces::coord_neighbor_fluid_fa7_gradU_1_
DoubleTab coord_neighbor_fluid_fa7_gradU_1_
Definition Post_Processing_Hydrodynamic_Forces.h:285
Post_Processing_Hydrodynamic_Forces::method_pressure_force_computation_
Method_pressure_force_computation method_pressure_force_computation_
Definition Post_Processing_Hydrodynamic_Forces.h:175
Post_Processing_Hydrodynamic_Forces::dVdz_P1_
DoubleTab dVdz_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:268
Post_Processing_Hydrodynamic_Forces::total_heat_transfer_
DoubleVect total_heat_transfer_
Definition Post_Processing_Hydrodynamic_Forces.h:234
Post_Processing_Hydrodynamic_Forces::dVdz_P2_
DoubleTab dVdz_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:278
Post_Processing_Hydrodynamic_Forces::pressure_force_fa7_
DoubleTab pressure_force_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:246
Post_Processing_Hydrodynamic_Forces::fill_gradU_P1
void fill_gradU_P1(int fa7, DoubleTab gradU_P1)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1699
Post_Processing_Hydrodynamic_Forces::prop_P2_fluid_compo_
DoubleTab prop_P2_fluid_compo_
Definition Post_Processing_Hydrodynamic_Forces.h:245
Post_Processing_Hydrodynamic_Forces::dVdx_P1_
DoubleTab dVdx_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:266
Post_Processing_Hydrodynamic_Forces::interpolation_distance_temperature_P1_
double interpolation_distance_temperature_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:224
Post_Processing_Hydrodynamic_Forces::interpolation_distance_gradU_P1_
double interpolation_distance_gradU_P1_
Definition Post_Processing_Hydrodynamic_Forces.h:228
Post_Processing_Hydrodynamic_Forces::interpolation_distance_gradU_P2_
double interpolation_distance_gradU_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:229
Post_Processing_Hydrodynamic_Forces::coord_neighbor_fluid_fa7_temp_1_
DoubleTab coord_neighbor_fluid_fa7_temp_1_
Definition Post_Processing_Hydrodynamic_Forces.h:287
Post_Processing_Hydrodynamic_Forces::proportion_fa7_ok_UP2_
DoubleTab proportion_fa7_ok_UP2_
Definition Post_Processing_Hydrodynamic_Forces.h:242
Post_Processing_Hydrodynamic_Forces::is_post_process_pressure_fa7_
bool is_post_process_pressure_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:217
Post_Processing_Hydrodynamic_Forces::resize_coord_neighbor_fluid_fa7
void resize_coord_neighbor_fluid_fa7(int nb_fa7)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1790
Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_gradU_elem_P1
int trilinear_interpolation_gradU_elem_P1(const DoubleTab &valeurs_champ, DoubleTab &coord, DoubleTab &resu)
Definition Post_Processing_Hydrodynamic_Forces.cpp:711
Post_Processing_Hydrodynamic_Forces::sigma_zx_fa7_
DoubleTab sigma_zx_fa7_
Definition Post_Processing_Hydrodynamic_Forces.h:259
Post_Processing_Hydrodynamic_Forces::dWdz_P2_
DoubleTab dWdz_P2_
Definition Post_Processing_Hydrodynamic_Forces.h:281
Post_Processing_Hydrodynamic_Forces::trilinear_interpolation_gradU_face
int trilinear_interpolation_gradU_face(const DoubleTab &valeurs_champ, DoubleTab &coord, DoubleTab &resu)
Definition Post_Processing_Hydrodynamic_Forces.cpp:501
Post_Processing_Hydrodynamic_Forces::list_elem_diph_
IntTab list_elem_diph_
Definition Post_Processing_Hydrodynamic_Forces.h:240
Post_Processing_Hydrodynamic_Forces::face_voisins_for_interp
int face_voisins_for_interp(int num_face, int i) const
Definition Post_Processing_Hydrodynamic_Forces.cpp:242
Post_Processing_Hydrodynamic_Forces::method_friction_force_computation_
Method_friction_force_computation method_friction_force_computation_
Definition Post_Processing_Hydrodynamic_Forces.h:177
Post_Processing_Hydrodynamic_Forces::Method_friction_force_computation
Method_friction_force_computation
Definition Post_Processing_Hydrodynamic_Forces.h:108
Post_Processing_Hydrodynamic_Forces::Method_friction_force_computation::TRILINEAR_LINEAR_PROJECTED_TENSOR
@ TRILINEAR_LINEAR_PROJECTED_TENSOR
Definition Post_Processing_Hydrodynamic_Forces.h:110
Post_Processing_Hydrodynamic_Forces::Method_friction_force_computation::TRILINEAR_LINEAR_COMPLET_TENSOR
@ TRILINEAR_LINEAR_COMPLET_TENSOR
Definition Post_Processing_Hydrodynamic_Forces.h:109
Post_Processing_Hydrodynamic_Forces::compute_viscosity_edges_sphere
double compute_viscosity_edges_sphere(int face1, int face2, int particle_id)
Definition Post_Processing_Hydrodynamic_Forces.cpp:1594
Process::mp_sum_for_each_item
static void mp_sum_for_each_item(TRUSTArray< _TYPE_ > &x, int n=-1)
Definition Process.cpp:193
Process::Journal
static Sortie & Journal(int message_level=0)
Renvoie un objet statique de type Sortie qui sert de journal d'evenements.
Definition Process.cpp:588
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Solid_Particle_base
: class Solid_Particle
Definition Solid_Particle_base.h:35
Solid_Particle_base::get_equivalent_radius
const double & get_equivalent_radius() const
Definition Solid_Particle_base.h:55
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133
Transport_Interfaces_FT_Disc
Definition Transport_Interfaces_FT_Disc.h:50
Transport_Interfaces_FT_Disc::maillage_interface_pour_post
virtual const Maillage_FT_Disc & maillage_interface_pour_post() const
Renvoie le maillage stocke specialement pour le postraitement (si on veut postraiter un etat intermed...
Definition Transport_Interfaces_FT_Disc.cpp:8980
Transport_Interfaces_FT_Disc::get_indicatrice
const Champ_base & get_indicatrice() override
getter champ variables_internes_->indicatrice_cache a partir de la position des interfaces.
Definition Transport_Interfaces_FT_Disc.cpp:1949
Transport_Interfaces_FT_Disc::maillage_interface
const Maillage_FT_Disc & maillage_interface() const
Definition Transport_Interfaces_FT_Disc.cpp:8336
Transport_Interfaces_FT_Disc::get_particles_velocity
const DoubleTab & get_particles_velocity() const
Definition Transport_Interfaces_FT_Disc.h:320
Transport_Interfaces_FT_Disc::get_particles_position
const DoubleTab & get_particles_position() const
Definition Transport_Interfaces_FT_Disc.h:319