16#include <Prepro_IBM_Ponderation.h>
17#include <Connectivite_som_elem.h>
25void Prepro_IBM_Ponderation::set_param(
Param& param)
const
40 Cout<<
"Weighting method = arimethic weight"<<finl;
42 Cout<<
"Weighting method = area weight"<<finl;
44 Cout<<
"Weighting method = inverse distance weight"<<finl;
46 Cout<<
"Weighting method = area and inverse distance weight"<<finl;
49 Cerr<<
"Prepro_IBM_ponderation : Type_de_ponderation : invalide argument = "<<pond_<<finl;
66 if (
verbose_) Cerr<<
"Prepro_IBM_Ponderation:: Computing arrays..."<<finl;
69 DoubleTab& normalArray = champ_normal_->valeurs();
71 DoubleTab& isNodeDirichletArray = isNodeDirichlet_->valeurs();
72 DoubleTab& rotationArray = champ_rotation_->valeurs();
73 int nbElemVol = rotationArray.
dimension(0);
78 if (maj_from_ext == 1)
80 isNodeDirichletArray = -1.;
82 for (
int elem = 0; elem <nbElemVol; elem++)
84 if (aireArray(elem) > 0.)
86 for (
int k = 0; k <dim_esp; k++) normalArray(elem,k) = rotationArray(elem,3*k+(dim_esp-1));
87 for (
int l = 0; l < le_dom_dis.
domaine().nb_som_elem(); l++)
89 isNodeDirichletArray(elems(elem,l)) = 1.0 ;
105 projectSolidPoints();
108 projectFluidPoints();
118void Prepro_IBM_Ponderation::projectSolidPoints()
120 const DoubleTab& normalArray = champ_normal_->valeurs();
121 const DoubleTab& aireArray =
champ_aire_->valeurs();
122 const DoubleTab& baryArray = champ_bary_->valeurs();
123 DoubleTab& solideArray = solid_points_->valeurs();
124 const DoubleTab& isNodeDirichletArray = isNodeDirichlet_->valeurs();
125 DoubleTab& solid_elemsArray = solid_elems_->valeurs();
127 solid_elemsArray = -2.;
129 Cout <<
"============================================================================"<<finl;
130 Cout<<
"prepro_IBM::projectSolidPoints - Weighting method = "<<pond_<<finl;
133 const Domaine& le_dom = le_dom_dis.
domaine();
134 const int nb_elem = le_dom.
nb_elem();
135 const int nb_som = le_dom.
nb_som();
137 const IntTab& connectDom = le_dom.
les_elems() ;
140 assert(dim_esp == 3);
142 DoubleTab x(dim_esp);
143 DoubleTrav sumContribProjSolidArray(nb_som);
144 sumContribProjSolidArray = 0.;
145 IntLists elemPerNeighbor(nb_som);
146 for (
int e=0; e<nb_elem; e++)
155 for (
int noeud=0; noeud<nb_som_elem; noeud++)
157 int sommet = connectDom(e,noeud);
159 x(0) = coordsDom3D(sommet,0) - baryArray(e,0);
163 x(1) = coordsDom3D(sommet,1) - baryArray(e,1);
167 x(2) = coordsDom3D(sommet,2) - baryArray(e,2);
175 double distance=0.0 , distance_normale=0.0;
176 for (
int d=0; d<dim_esp; d++)
178 distance_normale += x(d)*normalArray(e,d);
179 distance += x(d)*x(d);
181 distance = sqrt(distance);
195 pond = (distance>
eps_effec_)? (1.0/distance) : 1.0;
198 pond = (distance>
eps_effec_)? (aireArray(e)/distance) : aireArray(e);
206 for (
int d=0; d<dim_esp; d++) solideArray(sommet,d) += pond * (coordsDom3D(sommet,d)- distance_normale*normalArray(e,d));
207 sumContribProjSolidArray(sommet) += pond;
218 for (
int l=0; l<taille; l++)
220 int num_som_v= voisins[l];
221 if ( isNodeDirichletArray(num_som_v)<0 && !(elemPerNeighbor[num_som_v].contient(e)) )
225 elemPerNeighbor[num_som_v].add(e);
227 double dx = coordsDom3D(num_som_v, 0) - baryArray(e, 0);
228 double dy = coordsDom3D(num_som_v, 1) - baryArray(e, 1);
229 double dz = coordsDom3D(num_som_v, 2) - baryArray(e, 2);
231 double distance_normale_v = dx * normalArray(e, 0) + dy * normalArray(e, 1) + dz * normalArray(e, 2);
232 double d1 = sqrt(dx * dx + dy * dy + dz * dz);
241 pond_v = fabs(1.0 / d1);
244 pond_v = fabs(aireArray(e) / d1);
253 Cerr <<
"DIST_PROB : d1 ~ 0 for neighbor" << num_som_v <<
" of the element " << e << finl;
254 pond_v = (pond_ == 4) ? aireArray(e) : 1.0;
258 solideArray(num_som_v, 0) += pond_v * (coordsDom3D(num_som_v, 0) - distance_normale_v * normalArray(e, 0));
259 solideArray(num_som_v, 1) += pond_v * (coordsDom3D(num_som_v, 1) - distance_normale_v * normalArray(e, 1));
260 solideArray(num_som_v, 2) += pond_v * (coordsDom3D(num_som_v, 2) - distance_normale_v * normalArray(e, 2));
262 sumContribProjSolidArray(num_som_v) += pond_v;
269 for(
int sommet=0; sommet<nb_som; sommet++)
271 if (sumContribProjSolidArray(sommet))
273 for (
int d=0; d<dim_esp; d++) solideArray(sommet,d) /= sumContribProjSolidArray(sommet);
274 int interSoElem = le_dom.
chercher_elements(solideArray(sommet,0),solideArray(sommet,1),solideArray(sommet,2));
275 if (interSoElem == -1)
276 solid_elemsArray(sommet) = -1.;
278 solid_elemsArray(sommet) = float(interSoElem);
288void Prepro_IBM_Ponderation::projectFluidPoints()
291 DoubleTab& normalArray = champ_normal_->valeurs();
292 DoubleTab& fluide = fluid_points_->valeurs();
293 const DoubleTab& dirichlet = isNodeDirichlet_->valeurs();
294 const DoubleTab& hmax_node = h_max_node_->valeurs();
295 DoubleVect& elem_fluide = fluid_elems_->valeurs();;
298 const Domaine_dis_base& le_dom_dis = mon_pb_->domaine_dis();
299 const Domaine& le_dom = le_dom_dis.
domaine();
301 const IntTab& connect = le_dom.
les_elems();
303 const int nb_som = le_dom.
nb_som();
304 const int nb_sommets_tot = le_dom.
nb_som_tot();
305 const int nb_elem = le_dom.
nb_elem();
307 assert(dim_esp == 3);
309 DoubleTab normale_som(nb_som,dim_esp);
311 double x, y, z, xp, yp, zp, norme;
315 for (
int som = 0; som < nb_som; som++)
317 if (dirichlet(som) == 1.)
322 x = coordsDom3D(som, 0);
323 y = coordsDom3D(som, 1);
324 z = coordsDom3D(som, 2);
325 xp = solid_points_->valeurs()(som, 0);
326 yp = solid_points_->valeurs()(som, 1);
327 zp = solid_points_->valeurs()(som, 2);
328 norme = sqrt((x - xp)*(x - xp) + (y - yp)*(y - yp) + (z - zp)*(z - zp));
332 elem_fluide(som) = -3;
333 for (
int d = 0; d < dim_esp; d++)
335 fluide(som, d) = coordsDom3D(som, d);
336 normale_som(som, d) = 0.0;
341 for (
int d = 0; d < dim_esp; d++)
343 normale_som(som, d) = (coordsDom3D(som, d) - solid_points_->valeurs()(som, d)) / norme;
344 fluide(som, d) =
dimTab_(d) ? coordsDom3D(som, d) + fact * hmax_node(som) * normale_som(som, d) : coordsDom3D(som, d);
346 interFlElem = le_dom.
chercher_elements(fluide(som,0),fluide(som,1),fluide(som,2));
347 if (interFlElem == -1)
348 elem_fluide(som) = -3.;
351 elem_fluide(som) = float(interFlElem);
352 if (aire(interFlElem) > 0.0)
353 elem_fluide(som) = -1.;
357 for (
int j = 0; j < nb_som_elem; j++)
359 int s = connect(interFlElem, j);
360 if (dirichlet(s) == 1.0) flag =
true;
362 if (flag) elem_fluide(som) = -1.0;
368 if (
c_prepro_ > 0 && elem_fluide(som) == -1)
371 for (
int d = 0; d < dim_esp; d++)
372 fluide(som, d) =
dimTab_(d) ? coordsDom3D(som, d) + fact * hmax_node(som) * normale_som(som, d) : coordsDom3D(som, d);
373 interFlElem = le_dom.
chercher_elements(fluide(som,0),fluide(som,1),fluide(som,2));
374 if (interFlElem == -1)
375 elem_fluide(som) = -3.;
378 elem_fluide(som) = float(interFlElem);
379 if (aire(interFlElem) > 0.0)
380 elem_fluide(som) = -1.;
384 for (
int j = 0; j < nb_som_elem; j++)
386 int s = connect(interFlElem, j);
387 if (dirichlet(s) == 1.0) flag =
true;
389 if (flag) elem_fluide(som) = -1.0;
403 for (
int l=0; l<taille; l++)
405 int num_som_v= voisins[l];
406 if ( dirichlet(num_som_v)<0 )
408 x = coordsDom3D(num_som_v, 0);
409 y = coordsDom3D(num_som_v, 1);
410 z = coordsDom3D(num_som_v, 2);
411 xp = solid_points_->valeurs()(num_som_v, 0);
412 yp = solid_points_->valeurs()(num_som_v, 1);
413 zp = solid_points_->valeurs()(num_som_v, 2);
414 norme = sqrt((x - xp)*(x - xp) + (y - yp)*(y - yp) + (z - zp)*(z - zp));
418 elem_fluide(num_som_v) = -3;
419 for (
int d = 0; d < dim_esp; d++)
421 fluide(num_som_v, d) = coordsDom3D(som, d);
422 normale_som(num_som_v, d) = 0.0;
427 for (
int d = 0; d < dim_esp; d++)
429 normale_som(num_som_v, d) = (coordsDom3D(num_som_v, d) - solid_points_->valeurs()(num_som_v, d)) / norme;
430 fluide(num_som_v, d) =
dimTab_(d) ? coordsDom3D(num_som_v, d) + fact * hmax_node(num_som_v) * normale_som(num_som_v, d) : coordsDom3D(num_som_v, d);
432 interFlElem = le_dom.
chercher_elements(fluide(num_som_v,0),fluide(num_som_v,1),fluide(num_som_v,2));
433 if (interFlElem == -1)
434 elem_fluide(num_som_v) = -3.;
437 elem_fluide(num_som_v) = float(interFlElem);
438 if (aire(interFlElem) > 0.0)
439 elem_fluide(num_som_v) = -1.;
443 for (
int j = 0; j < nb_som_elem; j++)
445 int s = connect(interFlElem, j);
446 if (dirichlet(s) == 1.0) flag =
true;
448 if (flag) elem_fluide(num_som_v) = -1.0;
454 if (
c_prepro_ > 0 && elem_fluide(num_som_v) == -1)
457 for (
int d = 0; d < dim_esp; d++)
458 fluide(num_som_v, d) =
dimTab_(d) ? coordsDom3D(num_som_v, d) + fact * hmax_node(num_som_v) * normale_som(num_som_v, d) : coordsDom3D(num_som_v, d);
459 interFlElem = le_dom.
chercher_elements(fluide(num_som_v,0),fluide(num_som_v,1),fluide(num_som_v,2));
460 if (interFlElem == -1)
461 elem_fluide(num_som_v) = -3.;
464 elem_fluide(num_som_v) = float(interFlElem);
465 if (aire(interFlElem) > 0.0)
466 elem_fluide(num_som_v) = -1.;
470 for (
int j = 0; j < nb_som_elem; j++)
472 int s = connect(interFlElem, j);
473 if (dirichlet(s) == 1.0) flag =
true;
475 if (flag) elem_fluide(num_som_v) = -1.0;
488 Static_Int_Lists connectivite_som_elem;
489 construire_connectivite_som_elem(nb_sommets_tot, connect, connectivite_som_elem, 1 );
491 DoubleTrav normalWeighArray(nb_elem);
492 for (
int som = 0; som < nb_som; som++)
494 double n0 = normale_som(som,0);
495 double n1 = normale_som(som,1);
496 double n2 = normale_som(som,2);
497 norme = sqrt(n0*n0+n1*n1+n2*n2);
501 for (
int j = 0; j <sz ; j++)
503 int numElem = connectivite_som_elem(som,j);
507 normalArray(numElem,0) += n0;
508 normalArray(numElem,1) += n1;
509 normalArray(numElem,2) += n2;
510 normalWeighArray(numElem) += 1.0;
516 DoubleTrav t1Arr(nb_elem, dim_esp), t2Arr(nb_elem, dim_esp);
517 for (
int e=0; e<nb_elem; e++)
521 norme = sqrt(normalArray(e,0)*normalArray(e,0)+normalArray(e,1)*normalArray(e,1)+normalArray(e,2)*normalArray(e,2));
524 for (
int d = 0; d < dim_esp; d++) normalArray(e,d) /= norme;
int nb_som_elem() const
Renvoie le nombre de sommets des elements geometriques constituants le domaine.
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.
const DoubleTab_t & coord_sommets() const
int_t nb_som_tot() const
Renvoie le nombre total de sommets du domaine i.e. le nombre de sommets reels et virtuels sur le proc...
int_t nb_som() const
Renvoie le nombre de sommets du domaine.
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.
const Domaine & domaine() const
Class defining operators and methods for all reading operation in an input flow (file,...
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Helper class to factorize the readOn method of Objet_U classes.
void ajouter(const char *keyword, const int *value, Param::Nature nat=Param::OPTIONAL)
Register an integer parameter.
int lire_avec_accolades_depuis(Entree &is)
Parse the parameter block { ... } from is.
void compute_solid_fluid(int) override
void associer_pb(const Probleme_base &) override
void computeLocalFrame(const DoubleTab &, DoubleTab &, DoubleTab &, int)
void compute_h_max_elem()
void compute_NeighNode(int)
void set_param(Param &) const override
void computeMatRot(const DoubleTab &, DoubleTab &, DoubleTab &, int)
virtual void associer_pb(const Probleme_base &)
IntLists sommets_voisins_
classe Probleme_base C'est un Probleme_U qui n'est pas un couplage.
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Classe de base des flux de sortie.
int_t get_list_size(int_t i_liste) const
renvoie le nombre d'elements de la liste i
_SIZE_ dimension(int d) const
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")