TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
Champ_P1NC_implementation.cpp
1/****************************************************************************
2* Copyright (c) 2026, 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 <Champ_P1NC_implementation.h>
17#include <Dirichlet_paroi_defilante.h>
18#include <Dirichlet_paroi_fixe.h>
19#include <Modele_turbulence_hyd_base.h>
20#include <Champ_Fonc_P1NC.h>
21#include <Equation_base.h>
22#include <distances_VEF.h>
23#include <LecFicDiffuse.h>
24#include <Matrice_Bloc.h>
25#include <TRUSTLists.h>
26#include <Champ_P1NC.h>
27#include <Periodique.h>
28#include <Solv_GCP.h>
29#include <Domaine.h>
30#include <Debog.h>
31#include <SSOR.h>
32#include <kokkos++.h>
33#include <TRUSTArray_kokkos.tpp>
34
35static int methode_calcul_valeurs_sommets_static = 2 ;
36
37// -1 moyenne des faces
38// 1 moyenne des sommets sans min/max
39// 2 moyenne des sommets avec min/max (ancienne version)
40
41Champ_P1NC_implementation::Champ_P1NC_implementation()
42{
43 // Initialise l'attribut
44 filtrer_L2_deja_appele_=0;
45 // definition des types de solveurs par defaut (sinon, lecture dans solveur.bar)
46 solveur_L2.typer("Solv_GCP");
47 OWN_PTR(Precond_base) p;
48 p.typer("SSOR");
49 ref_cast(Solv_GCP,solveur_L2.valeur()).set_precond(p);
50 solveur_L2.nommer("solveur_L2");
51 solveur_H1.typer("Solv_GCP");
52 ref_cast(Solv_GCP,solveur_H1.valeur()).set_precond(p);
53 solveur_H1.nommer("solveur_H1");
54}
55
56static const double coeff_penalisation = 1e9;
57
58/*! @brief Projection du champ P1NC "cha" sur l'espace des champs P1.
59 *
60 * Le resultat est stocke dans cha.ch_som() et renvoye (valeur de retour)
61 * Voir note technique 2006/010 de Patrick Quemere "Nouvelle approche VEF pour la LES...".
62 * Voir aussi "What is in TRUST" : Scales_Separation_P1_P1NC.
63 *
64 */
65DoubleTab& valeur_P1_L2(Champ_P1NC& cha, const Domaine& dom)
66{
67 int nb_elem_tot = dom.nb_elem_tot();
68 int nb_som = dom.nb_som_tot();
69 int nb_som_elem = dom.nb_som_elem();
70 const DoubleTab& ch = cha.valeurs();
71
72 const Domaine_VEF& domaine_VEF = cha.domaine_vef();
73 const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF, cha.equation().domaine_Cl_dis());
74 const IntTab& elem_som = dom.les_elems();
75 double mijK, miiK;
76 DoubleTab& vit_som=cha.ch_som();
77
78 int n_bord, nb_comp=vit_som.line_size();
79 DoubleTab secmem(cha.ch_som());
80
81 secmem=0.;
82 // vit_som=0.; // PQ : 06/10/04 : dans la resolution du systeme, "solution" repart de vit_som=cha.ch_som()
83 // pour s'affranchir de recalculer u_bar lors d'un double appel a filtrer_L2
84
85 if(Objet_U::dimension==2)
86 {
87 mijK=1./12.;
88 miiK= 1./6.;
89 }
90 else
91 {
92 mijK=1./18.;
93 miiK= 1./12.;
94 }
95
96 // dans le cas d'un domaine nul on doit effectuer le dimensionnement
97 // mais une seule fois
98 double non_prepare=0;
99 if (cha.MatP1NC2P1_L2.nb_lignes() <2)
100 non_prepare=1;
101 non_prepare=Process::mp_min(non_prepare);
102 if (non_prepare==1)
103 {
104 // Initialisation du solveur
105 SolveurSys& solv=ref_cast(SolveurSys, (cha.solveur_L2));
106 LecFicDiffuse fic;
107 fic.ouvrir("solveur.bar");
108 if(fic.good())
109 {
110 fic >> solv;
111 Cerr << Process::me() << "solveur.bar ouvert pour lire: " << solv << finl;
112 }
113 solv->reinit();
114
115 Cerr<<"On remplit la matrice cha.MatP1NC2P1_L2 pour le champ "<< cha.le_nom() <<finl;
116 Matrice_Morse_Sym& mat=ref_cast(Matrice_Morse_Sym, (cha.MatP1NC2P1_L2));
117 int rang;
118 int nnz=0;
119 IntLists voisins(nb_som);
120 DoubleLists coeffs(nb_som);
121 DoubleVect diag(nb_som);
122 int elem;
123 for (elem=0; elem<nb_elem_tot; elem++)
124 {
125 double coeff_ij=mijK*domaine_VEF.volumes(elem);
126 double coeff_ii=miiK*domaine_VEF.volumes(elem);
127 for (int isom=0; isom<nb_som_elem; isom++)
128 {
129 int som=elem_som(elem,isom);
130 int ii=dom.get_renum_som_perio(som);
131 diag[ii]+=coeff_ii;
132 for (int jsom=isom+1; jsom<nb_som_elem; jsom++)
133 {
134 int asom=elem_som(elem,jsom);
135 int j=dom.get_renum_som_perio(asom);
136 int i=ii;
137 if(i>j)
138 {
139 int k=j;
140 j=i;
141 i=k;
142 }
143 rang=voisins[i].rang(j);
144 if(rang==-1)
145 {
146 voisins[i].add(j);
147 coeffs[i].add(coeff_ij);
148 nnz++;
149 }
150 else
151 {
152 coeffs[i][rang]+=coeff_ij;
153 }
154 }
155 }
156 }
157
158 // PQ : pour le remplissage des coefficients de la matrice specifiques au CL Dirichlet,
159 // on commence par traiter les conditions Dirichlet, PUIS Dirichlet_homogene pour s'assurer
160 // de l'"adherence" des points sommets appartenant aux deux types de CL (aretes de bords)
161
162 IntVect test_cl_imposee(nb_som);
163 test_cl_imposee = 0;
164
165 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
166 {
167 //for n_bord
168 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
169 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
170 int face;
171 int num1 = 0;
172 int num2 = le_bord.nb_faces_tot();
173 // if ( sub_type(Dirichlet,la_cl.valeur()))
174 if (!sub_type(Periodique,la_cl.valeur()))
175 {
176 for (int ind_face=num1; ind_face<num2; ind_face++)
177 {
178 face = le_bord.num_face(ind_face);
179 for(int isom=0; isom<Objet_U::dimension; isom++)
180 {
181 int som=domaine_VEF.face_sommets(face,isom);
182 som=dom.get_renum_som_perio(som);
183 if (test_cl_imposee(som) == 0)
184 diag[som] = 0;
185 diag[som] += coeff_penalisation;
186 test_cl_imposee(som) += 1;
187 }
188 }
189 }
190 }//for n_bord
191
192 test_cl_imposee = 0;
193
194 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
195 {
196 //for n_bord
197 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
198 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
199 int face;
200 int num1 = 0;
201 int num2 = le_bord.nb_faces_tot();
202
203 if ( sub_type(Dirichlet_homogene,la_cl.valeur()) )
204 {
205 for (int ind_face=num1; ind_face<num2; ind_face++)
206 {
207 face = le_bord.num_face(ind_face);
208 for(int isom=0; isom<Objet_U::dimension; isom++)
209 {
210 int som=domaine_VEF.face_sommets(face,isom);
211 som=dom.get_renum_som_perio(som);
212 if (test_cl_imposee(som) == 0)
213 diag[som] = 0;
214 diag[som] += coeff_penalisation;
215 test_cl_imposee(som) += 1;
216 }
217 }
218 }
219 }//for n_bord
220
221 {
222 for (int i=0; i<nb_som; i++)
223 if(diag(i)==0)
224 diag(i)=1.;
225 }
226
227 // coeffs = 0.;
228 // diag = 10.;
229 mat.dimensionner(nb_som, nnz+nb_som) ;
230 mat.remplir(voisins, coeffs, diag) ;
231 mat.compacte() ;
232 mat.set_est_definie(1);
233
234 // Modif CD 15/04/03 pour paralleliser le filtrage
235 // On copie la matrice remplie ci-dessus dans une matrice bloc.
236 // Cette modification permet d'utiliser le preconditionnement SSOR
237
238 //On fait aussi la transformation en matrice bloc en sequentiel car sinon
239 //probleme avec NP
240
241 cha.dimensionner_Mat_Bloc_Morse_Sym(cha.MatP1NC2P1_L2_Parallele);
242 cha.Mat_Morse_to_Mat_Bloc(cha.MatP1NC2P1_L2_Parallele);
243
244 // Fin Modif CD 15/04/03
245
246 //mat.imprimer_formatte(Cerr);
247 // Cerr<<"domaine.nb_som_tot() = "<<dom.nb_som_tot()<<", domaine.nb_som() = "<<dom.nb_som()<<finl;
248 // Cerr << "Matrice de filtrage OK" << finl;
249 // Cout << "Matrice de masse P1 : " << finl;
250 // mat.imprimer(Cout) << finl;
251 }
252
253 // Debog::verifier(" verif champ " , ch);
254
255 double coeff;
256 if(Objet_U::dimension==2)
257 coeff=1./6.;
258 else
259 coeff=1./12.;
260
261 int nb_faces=domaine_VEF.nb_faces();
262 int nb_faces_tot=domaine_VEF.nb_faces_tot();
263 ArrOfInt virtfait(nb_faces_tot-nb_faces);
264 int nb_som_face=domaine_VEF.nb_som_face();
265 const IntTab& face_sommets=domaine_VEF.face_sommets();
266 const IntTab& face_voisins=domaine_VEF.face_voisins();
267 int face;
268 int premiere_face_int = domaine_VEF.premiere_face_int();
269
270 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
271 {
272 //for n_bord
273 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
274 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
275 int nb_faces_bord = le_bord.nb_faces();
276 int num1 = 0;
277 int num2 = le_bord.nb_faces_tot();
278
279 if (sub_type(Periodique,la_cl.valeur()))
280 {
281 // *0.5 car les faces perio sont parcourues deux fois
282 for (int ind_face=num1; ind_face<num2; ind_face++)
283 {
284 face = le_bord.num_face(ind_face);
285 if (ind_face-nb_faces_bord>=0) virtfait[face-nb_faces]=1;
286 int elem1=face_voisins(face,0);
287 int elem2=face_voisins(face,1);
288 double volume=domaine_VEF.volumes(elem1);
289 if(elem2!=-1)
290 volume+=domaine_VEF.volumes(elem2);
291 for(int isom=0; isom<nb_som_face; isom++)
292 {
293 int som=face_sommets(face,isom);
294 int som1=dom.get_renum_som_perio(som);
295 for(int k=0; k<nb_comp; k++)
296 secmem(som1,k)+=0.5*coeff*volume*ch(face,k);
297 }
298 }
299 }
300 else
301 for (int ind_face=num1; ind_face<num2; ind_face++)
302 {
303 face = le_bord.num_face(ind_face);
304 if (ind_face-nb_faces_bord>=0) virtfait[face-nb_faces]=1;
305 int elem1=face_voisins(face,0);
306 int elem2=face_voisins(face,1);
307 double volume=domaine_VEF.volumes(elem1);
308 if(elem2!=-1)
309 volume+=domaine_VEF.volumes(elem2);
310 for(int isom=0; isom<nb_som_face; isom++)
311 {
312 int som=face_sommets(face,isom);
313 int som1=dom.get_renum_som_perio(som);
314 for(int k=0; k<nb_comp; k++)
315 secmem(som1,k)+=coeff*volume*ch(face,k);
316 }
317 }
318 }
319
320 for(face=premiere_face_int; face<nb_faces_tot; face++)
321 {
322 int indvirtface=face-nb_faces;
323 // on regarde si la face virtuelle n'a pas deja ete traitee
324 if (indvirtface>=0)
325 if (virtfait[indvirtface]==1) continue;
326
327 int elem1=face_voisins(face,0);
328 int elem2=face_voisins(face,1);
329 double volume=domaine_VEF.volumes(elem1);
330 if(elem2!=-1)
331 volume+=domaine_VEF.volumes(elem2);
332 for(int isom=0; isom<nb_som_face; isom++)
333 {
334 int som=face_sommets(face,isom);
335 int som1=dom.get_renum_som_perio(som);
336 for(int k=0; k<nb_comp; k++)
337 secmem(som1,k)+=coeff*volume*ch(face,k);
338 }
339 }
340 // ************************************************************************************************
341 // ************************************************************************************************
342 //
343 // Pour imposer u_som = u donne par CL_Dirichlet
344 ////
345 // ************************************************************************************************
346 // ************************************************************************************************
347
348 // PQ : on commence par traiter les conditions Dirichlet, PUIS Dirichlet_homogene pour s'assurer
349 // de l'"adherence" des points sommets appartenant aux deux types de CL (aretes de bords)
350
351 IntVect test_cl_imposee(nb_som);
352
353 test_cl_imposee = 0;
354
355 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
356 {
357 //for n_bord
358 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
359 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
360 int num1 = 0;
361 int num2 = le_bord.nb_faces_tot();
362
363 // if ( sub_type(Dirichlet,la_cl.valeur()))
364 if (!sub_type(Periodique,la_cl.valeur()))
365 {
366 for (int ind_face=num1; ind_face<num2; ind_face++)
367 {
368 face = le_bord.num_face(ind_face);
369 for(int isom=0; isom<nb_som_face; isom++)
370 {
371 int som=face_sommets(face,isom);
372 int som1=dom.get_renum_som_perio(som);
373 for(int k=0; k<nb_comp; k++)
374 {
375 if (test_cl_imposee(som1)==0)
376 secmem(som1,k) = 0;
377
378 secmem(som1,k) += coeff_penalisation * ch(face,k);
379 }
380 test_cl_imposee(som1) += 1;
381 }
382 }
383 }
384 }//for n_bord
385
386 test_cl_imposee = 0;
387
388 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
389 {
390 //for n_bord
391 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
392 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
393 int num1 = 0;
394 int num2 = le_bord.nb_faces_tot();
395
396 if (sub_type(Dirichlet_homogene,la_cl.valeur()))
397 {
398 for (int ind_face=num1; ind_face<num2; ind_face++)
399 {
400 face = le_bord.num_face(ind_face);
401 for(int isom=0; isom<nb_som_face; isom++)
402 {
403 int som=face_sommets(face,isom);
404 int som1=dom.get_renum_som_perio(som);
405 for(int k=0; k<nb_comp; k++)
406 {
407 if (test_cl_imposee(som1)==0)
408 secmem(som1,k) = 0;
409
410 secmem(som1,k) += coeff_penalisation * ch(face,k);
411 }
412 test_cl_imposee(som1) += 1;
413 }
414 }
415 }
416 }//for n_bord
417
418 // ************************************************************************************************
419 // ************************************************************************************************
420 secmem.echange_espace_virtuel();
421 // Debog::verifier(" verif secmem_ech " , secmem);
422
423 SolveurSys& solv=ref_cast(SolveurSys, cha.solveur_L2);
424 DoubleVect secmemk(cha.ch_som_vect());
425 DoubleVect solution(cha.ch_som_vect());
426
427 if(nb_comp==1)
428 {
429 solution=vit_som;
430 secmem.echange_espace_virtuel();
431 Debog::verifier(" verif secmem " , secmem);
432
433
434 //ON UTILISE LA MEME MATRICE EN SEQUENTIEL ET PARALLELE
435 //CE QUI PERMET D UTILISER UNE MATRICE BLOC QUAND ON APPELLE LE SOLVEUR NP
436 solv.resoudre_systeme(cha.MatP1NC2P1_L2_Parallele.valeur(), secmem, solution);
437
438 for(int som=0; som<nb_som; som++)
439 vit_som(som)=solution(dom.get_renum_som_perio(som));
440 }
441 else
442 {
443 for(int k=0; k<nb_comp; k++)
444 {
445 int som;
446 for(som=0; som<nb_som; som++)
447 {
448 solution(som)=vit_som(som,k);
449 secmemk(som)=secmem(som,k);
450 }
451
452 solution.echange_espace_virtuel();
453 secmemk.echange_espace_virtuel();
454
455 //ON UTILISE LA MEME MATRICE EN SEQUENTIEL ET PARALLELE
456 //CE QUI PERMET D UTILISER UNE MATRICE BLOC QUAND ON APPELLE LE SOLVEUR NP
457 solv.resoudre_systeme(cha.MatP1NC2P1_L2_Parallele.valeur(), secmemk, solution);
458
459 for(som=0; som<nb_som; som++)
460 vit_som(som,k)=solution(dom.get_renum_som_perio(som));
461 }
462 }
463
464 // Debog::verifier(" verif vit_som 0 " , vit_som);
465 vit_som.echange_espace_virtuel();
466 // Debog::verifier(" verif vit_som 1 " , vit_som);
467
468 // {
469 // for(int k=0; k<nb_comp; k++)
470 // for(int som=0; som<nb_som; som++)
471 // Cerr << k << " x " << coord(som,0)<< " y " << coord(som,1) << " " << vit_som(som,k) << finl;
472 // double psc=vit_som*vit_som;
473 // psc = Process::mp_sum(psc);
474 // Debog::verifier("norme : ", psc);
475 // Cerr << "############ psc = " << psc << finl;
476 // }
477 Debog::verifier(" verif vit_som " , vit_som);
478 return vit_som;
479}
480
481DoubleTab& valeur_P1_L2(Champ_Fonc_P1NC& cha, const Domaine& dom)
482{
483 int nb_elem_tot = dom.nb_elem_tot();
484 int nb_som = dom.nb_som_tot();
485 int nb_som_elem = dom.nb_som_elem();
486 const DoubleTab& ch = cha.valeurs();
487 const Domaine_VEF& domaine_VEF = cha.domaine_vef();
488 const IntTab& elem_som = dom.les_elems();
489 double mijK, miiK;
490 DoubleTab& vit_som=cha.ch_som();
491
492 int nb_comp=vit_som.line_size();
493 DoubleTab secmem(cha.ch_som());
494 secmem=0.;
495
496 if(Objet_U::dimension==2)
497 {
498 mijK=1./12.;
499 miiK= 1./6.;
500 }
501 else
502 {
503 mijK=1./18.;
504 miiK= 1./12.;
505 }
506
507 // dans le cas d'un domaine nul on doit effectuer le dimensionnement
508 // mais une seule fois
509 double non_prepare=0;
510 if (cha.MatP1NC2P1_L2.nb_lignes() <2)
511 non_prepare=1;
512 non_prepare=Process::mp_min(non_prepare);
513 if (non_prepare==1)
514 {
515 // Initialisation du solveur
516 SolveurSys& solv=ref_cast(SolveurSys, (cha.solveur_L2));
517 LecFicDiffuse fic;
518 fic.ouvrir("solveur.bar");
519 if(fic.good())
520 {
521 fic >> solv;
522 Cerr << Process::me() << "solveur.bar ouvert pour lire: " << solv << finl;
523 }
524 solv->reinit();
525
526 Cerr<<"On remplit la matrice cha.MatP1NC2P1_L2 pour le champ "<< cha.le_nom() <<finl;
527 Matrice_Morse_Sym& mat=ref_cast(Matrice_Morse_Sym, (cha.MatP1NC2P1_L2));
528 int rang;
529 int nnz=0;
530 IntLists voisins(nb_som);
531 DoubleLists coeffs(nb_som);
532 DoubleVect diag(nb_som);
533 int elem;
534 for (elem=0; elem<nb_elem_tot; elem++)
535 {
536 double coeff_ij=mijK*domaine_VEF.volumes(elem);
537 double coeff_ii=miiK*domaine_VEF.volumes(elem);
538 for (int isom=0; isom<nb_som_elem; isom++)
539 {
540 int som=elem_som(elem,isom);
541 int ii=dom.get_renum_som_perio(som);
542 diag[ii]+=coeff_ii;
543 for (int jsom=isom+1; jsom<nb_som_elem; jsom++)
544 {
545 int asom=elem_som(elem,jsom);
546 int j=dom.get_renum_som_perio(asom);
547 int i=ii;
548 if(i>j)
549 {
550 int k=j;
551 j=i;
552 i=k;
553 }
554 rang=voisins[i].rang(j);
555 if(rang==-1)
556 {
557 voisins[i].add(j);
558 coeffs[i].add(coeff_ij);
559 nnz++;
560 }
561 else
562 {
563 coeffs[i][rang]+=coeff_ij;
564 }
565 }
566 }
567 }
568
569
570 {
571 for (int i=0; i<nb_som; i++)
572 if(diag(i)==0)
573 diag(i)=1.;
574 }
575
576 mat.dimensionner(nb_som, nnz+nb_som) ;
577 mat.remplir(voisins, coeffs, diag) ;
578 mat.compacte() ;
579 mat.set_est_definie(1);
580
581 // Modif CD 15/04/03 pour paralleliser le filtrage
582 // On copie la matrice remplie ci-dessus dans une matrice bloc.
583 // Cette modification permet d'utiliser le preconditionnement SSOR
584
585 //On fait aussi la transformation en matrice bloc en sequentiel car sinon
586 //probleme avec NP
587
588 cha.dimensionner_Mat_Bloc_Morse_Sym(cha.MatP1NC2P1_L2_Parallele);
589 cha.Mat_Morse_to_Mat_Bloc(cha.MatP1NC2P1_L2_Parallele);
590
591 // Fin Modif CD 15/04/03
592
593 //mat.imprimer_formatte(Cerr);
594 // Cerr<<"domaine.nb_som_tot() = "<<dom.nb_som_tot()<<", domaine.nb_som() = "<<dom.nb_som()<<finl;
595 // Cerr << "Matrice de filtrage OK" << finl;
596 // Cout << "Matrice de masse P1 : " << finl;
597 // mat.imprimer(Cout) << finl;
598 }
599
600 // Debog::verifier(" verif champ " , ch);
601
602 double coeff;
603 if(Objet_U::dimension==2)
604 coeff=1./6.;
605 else
606 coeff=1./12.;
607
608 int nb_faces=domaine_VEF.nb_faces();
609 int nb_faces_tot=domaine_VEF.nb_faces_tot();
610 ArrOfInt virtfait(nb_faces_tot-nb_faces);
611 int nb_som_face=domaine_VEF.nb_som_face();
612 const IntTab& face_sommets=domaine_VEF.face_sommets();
613 const IntTab& face_voisins=domaine_VEF.face_voisins();
614 int face;
615
616 for(face=0; face<nb_faces_tot; face++)
617 {
618 int indvirtface=face-nb_faces;
619 // on regarde si la face virtuelle n'a pas deja ete traitee
620 if (indvirtface>=0)
621 if (virtfait[indvirtface]==1) continue;
622
623 int elem1=face_voisins(face,0);
624 int elem2=face_voisins(face,1);
625 double volume=domaine_VEF.volumes(elem1);
626 if(elem2!=-1)
627 volume+=domaine_VEF.volumes(elem2);
628 for(int isom=0; isom<nb_som_face; isom++)
629 {
630 int som=face_sommets(face,isom);
631 int som1=dom.get_renum_som_perio(som);
632 //Pour un sommet periodique on ne retient qu une seule contribution
633 //pour deux faces en vis a vis
634 if (som1==som)
635 for(int k=0; k<nb_comp; k++)
636 secmem(som1,k)+=coeff*volume*ch(face,k);
637 }
638 }
639
640 secmem.echange_espace_virtuel();
641 // Debog::verifier(" verif secmem_ech " , secmem);
642
643 SolveurSys& solv=ref_cast(SolveurSys, cha.solveur_L2);
644 DoubleVect secmemk(cha.ch_som_vect());
645 DoubleVect solution(cha.ch_som_vect());
646
647 if(nb_comp==1)
648 {
649 solution=vit_som;
650 secmem.echange_espace_virtuel();
651 Debog::verifier(" verif secmem " , secmem);
652
653
654 //ON UTILISE LA MEME MATRICE EN SEQUENTIEL ET PARALLELE
655 //CE QUI PERMET D UTILISER UNE MATRICE BLOC QUAND ON APPELLE LE SOLVEUR NP
656 solv.resoudre_systeme(cha.MatP1NC2P1_L2_Parallele.valeur(), secmem, solution);
657
658 for(int som=0; som<nb_som; som++)
659 vit_som(som)=solution(dom.get_renum_som_perio(som));
660 }
661 else
662 {
663 for(int k=0; k<nb_comp; k++)
664 {
665 int som;
666 for(som=0; som<nb_som; som++)
667 {
668 solution(som)=vit_som(som,k);
669 secmemk(som)=secmem(som,k);
670 }
671
672 solution.echange_espace_virtuel();
673 secmemk.echange_espace_virtuel();
674
675 //ON UTILISE LA MEME MATRICE EN SEQUENTIEL ET PARALLELE
676 //CE QUI PERMET D UTILISER UNE MATRICE BLOC QUAND ON APPELLE LE SOLVEUR NP
677 solv.resoudre_systeme(cha.MatP1NC2P1_L2_Parallele.valeur(), secmemk, solution);
678
679 for(som=0; som<nb_som; som++)
680 vit_som(som,k)=solution(dom.get_renum_som_perio(som));
681 }
682 }
683
684 // Debog::verifier(" verif vit_som 0 " , vit_som);
685 vit_som.echange_espace_virtuel();
686 // Debog::verifier(" verif vit_som 1 " , vit_som);
687
688 // {
689 // for(int k=0; k<nb_comp; k++)
690 // for(int som=0; som<nb_som; som++)
691 // Cerr << k << " x " << coord(som,0)<< " y " << coord(som,1) << " " << vit_som(som,k) << finl;
692 // double psc=vit_som*vit_som;
693 // psc = Process::mp_sum(psc);
694 // Debog::verifier("norme : ", psc);
695 // Cerr << "############ psc = " << psc << finl;
696 // }
697 Debog::verifier(" verif vit_som " , vit_som);
698 return vit_som;
699}
700
701DoubleTab&
702valeur_P1_H1(const Champ_P1NC& cha,
703 const Domaine& dom,
704 DoubleTab& ch_som)
705{
706 const Domaine& domaine = dom;
707 int nb_elem_tot = domaine.nb_elem_tot();
708 int nb_som = domaine.nb_som();
709 int nb_som_elem = domaine.nb_som_elem();
710 const DoubleTab& ch = cha.valeurs();
711 const Domaine_VEF& domaine_VEF = cha.domaine_vef();
712 const DoubleTab& normales=domaine_VEF.face_normales();
713 const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF, cha.equation().domaine_Cl_dis());
714 const IntTab& elem_som = domaine.les_elems();
715 const IntTab& elem_faces=domaine_VEF.elem_faces();
716 double mijK;
717 int nb_comp=ch_som.dimension(1);
718 int dimension=Objet_U::dimension;
719
720 if(dimension==2)
721 {
722 mijK=1./4.;
723 }
724 else
725 {
726 mijK=1./9.;
727 }
728
729 int elem;
730 int n_bord;
731 static double diag_ref;
732 double non_prepare=0;
733 if ( cha.MatP1NC2P1_H1.nb_lignes() <2 )
734 non_prepare=1;
735 non_prepare=Process::mp_min(non_prepare);
736 if (non_prepare==1)
737 {
738 EFichier fic;
739 if(fic.ouvrir("solveur.bar"))
740 {
741 SolveurSys& solv=ref_cast_non_const(SolveurSys,cha.solveur_H1);
742 fic >> solv;
743 }
744 Matrice_Morse_Sym& mat=ref_cast_non_const(Matrice_Morse_Sym,cha.MatP1NC2P1_H1);
745 int rang;
746 int nnz=0;
747 IntLists voisins(nb_som);
748 DoubleLists coeffs(nb_som);
749 DoubleVect diag(nb_som);
750 for (elem=0; elem<nb_elem_tot; elem++)
751 {
752 double volume=domaine_VEF.volumes(elem);
753 for (int isom=0; isom<nb_som_elem; isom++)
754 {
755 int facei=elem_faces(elem,isom);
756 int i=dom.get_renum_som_perio(elem_som(elem,isom));
757 for (int jsom=isom+1; jsom<nb_som_elem; jsom++)
758 {
759 int facej=elem_faces(elem,jsom);
760 int j=dom.get_renum_som_perio(elem_som(elem,jsom));
761 double coeffij=0.0;
762 for(int k=0; k<dimension; k++)
763 coeffij+=normales(facei,k)*normales(facej,k);
764
765 coeffij*=domaine_VEF.oriente_normale(facei,elem)*
766 domaine_VEF.oriente_normale(facej,elem);
767 coeffij*=mijK/volume;
768 int ii=i, jj=j;
769 if(ii>jj)
770 {
771 int tmp=ii;
772 ii=jj;
773 jj=tmp;
774 }
775 rang=voisins[ii].rang(jj);
776 if(rang==-1)
777 {
778 voisins[ii].add(jj);
779 coeffs[ii].add(coeffij);
780 nnz++;
781 }
782 else
783 {
784 coeffs[ii][rang]+=coeffij;
785 }
786 diag[ii]-=coeffij;
787 diag[jj]-=coeffij;
788 }
789 }
790 }
791 diag_ref=diag(0);
792 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
793 {
794 //for n_bord
795 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
796 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
797 int num1 = le_bord.num_premiere_face();
798 int num2 = num1 + le_bord.nb_faces();
799
800 for (int face=num1; face<num2; face++)
801 {
802 elem=domaine_VEF.face_voisins(face,0);
803 for(int isom=0; isom<dimension; isom++)
804 {
805 int som=domaine_VEF.face_sommets(face,isom);
806 som=dom.get_renum_som_perio(som);
807 diag[som]=1.e1*diag_ref;
808 }
809 }
810 }
811
812 mat.dimensionner(nb_som, nnz+nb_som) ;
813 mat.remplir(voisins, coeffs, diag) ;
814 mat.compacte() ;
815 mat.set_est_definie(1);
816 // Cerr << "Matrice de filtrage OK" << finl;
817 // Cerr << "Matrice de Laplacien P1 : " << finl;
818 // mat.imprimer(Cerr) << finl;
819 }
820
821
822 double coeff=-1./dimension;
823 int k;
824 for(int comp=0; comp<nb_comp; comp++)
825 {
826 //for comp
827 DoubleVect secmem(nb_som);
828 DoubleVect solution(nb_som);
829 for (elem=0; elem<nb_elem_tot; elem++)
830 {
831 //for elem
832 double volume=domaine_VEF.volumes(elem);
833 for (int isom=0; isom<nb_som_elem; isom++)
834 {
835 int facei=elem_faces(elem,isom);
836 int i=dom.get_renum_som_perio(elem_som(elem,isom));
837 for (int jsom=0; jsom<nb_som_elem; jsom++)
838 {
839 int facej=elem_faces(elem,jsom);
840 double coeffij=0;
841 for(k=0; k<dimension; k++)
842 coeffij+=normales(facei,k)*normales(facej,k);
843 coeffij*=domaine_VEF.oriente_normale(facei,elem)*
844 domaine_VEF.oriente_normale(facej,elem);
845 coeffij*=coeff/volume;
846 secmem(i)+=coeffij*ch(facej,comp);
847 }
848 }
849 }
850 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
851 {
852 //for n_bord
853 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
854 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
855 int num1 = le_bord.num_premiere_face();
856 int num2 = num1 + le_bord.nb_faces();
857
858 for (int face=num1; face<num2; face++)
859 {
860 elem=domaine_VEF.face_voisins(face,0);
861 for(int isom=0; isom<dimension; isom++)
862 {
863 int som=domaine_VEF.face_sommets(face,isom);
864 som=dom.get_renum_som_perio(som);
865 secmem(som)=1.e1*diag_ref* cha.valeur_a_sommet_compo(som, elem, comp);
866 }
867 }
868 }
869
870 SolveurSys& solv=ref_cast_non_const(SolveurSys,cha.solveur_H1);
871 solv.resoudre_systeme(cha.MatP1NC2P1_H1, secmem, solution);
872
873 for(int som=0; som<nb_som; som++)
874 ch_som(som,comp)=solution(dom.get_renum_som_perio(som));
875
876 }
877
878 return ch_som;
879}
880
881void Champ_P1NC_implementation::filtrer_H1(DoubleTab& valeurs) const
882{
883 //Cerr << "Champ_P1NC_implementation::filtrer" << finl;
884 const Domaine_VEF& domaine_VEF = domaine_vef();
885 static int test=1;
886 if(!test)
887 {
888 test=1;
889 DoubleTab xv=domaine_VEF.xv();
890 Cerr << "TEST filtrer_H1 solution = " << mp_norme_vect(xv) << finl;
891 filtrer_H1(xv);
892 xv-=domaine_VEF.xv();
893 Cerr << "TEST filtrer_H1 erreur = " << mp_norme_vect(xv) << finl;
894 Cerr << "err = " << xv
895 << "domaine_VEF.xv() " << domaine_VEF.xv() << finl;
896 xv=1.;
897 filtrer_H1(xv);
898 xv-=1.;
899 Cerr << "erreur sur les constantes : " << mp_norme_vect(xv) << finl;
900 }
901 int nb_som= domaine_VEF.domaine().nb_som_tot();
902 int nb_faces=valeurs.dimension(0);
903 const Champ_P1NC& cha_const =ref_cast(Champ_P1NC,le_champ());
904 Champ_P1NC& cha=ref_cast_non_const(Champ_P1NC,cha_const);
905 DoubleTab valeurs_fac(cha.valeurs());
906 cha.valeurs()=valeurs;
907 const int nb_comp=valeurs.line_size();
908 DoubleTab vit_som(nb_som, nb_comp);
909 valeur_P1_H1(cha, cha.domaine(), vit_som);
910 cha.valeurs()=valeurs_fac;
911 double coeff= 1./Objet_U::dimension;
912
913 for(int face=0; face<nb_faces; face++)
914 {
915 int s1=domaine_VEF.face_sommets(face,0);
916 int s2=domaine_VEF.face_sommets(face,1);
917 int s3=-1;
918 if(Objet_U::dimension==3)
919 s3=domaine_VEF.face_sommets(face,2);
920 for(int comp=0; comp<nb_comp; comp++)
921 {
922 double somme=0;
923 somme+=(vit_som(s1,comp)+vit_som(s2,comp));
924 if(Objet_U::dimension==3)
925 somme+=vit_som(s3,comp);
926 valeurs(face, comp)=coeff*somme;
927 }
928 }
929}
930
931void Champ_P1NC_implementation::filtrer_L2(DoubleTab& valeurs) const
932{
933 const Domaine_VEF& domaine_VEF = domaine_vef();
934 const Domaine& dom=domaine_VEF.domaine();
935 double coeff= 1./Objet_U::dimension;
936 // int nb_som= domaine_VEF.domaine().nb_som_tot();
937 int nb_faces=valeurs.dimension_tot(0);
938 Champ_base& cha = ref_cast_non_const(Champ_base,le_champ());
939 if (!sub_type(Champ_P1NC,cha) && (!sub_type(Champ_Fonc_P1NC,cha)))
940 {
941 Cerr<<"on ne doit pas filtre_L2 des champs qui ne sont pas de type P1NC"<<finl;
942 Process::exit();
943 }
944
945 const int nb_comp=valeurs.line_size();
946
947 if (sub_type(Champ_P1NC,cha))
948 {
949 Champ_P1NC& ch_inc_P1NC = ref_cast(Champ_P1NC,cha);
950
951 if (ch_inc_P1NC.equation().le_nom()!="gradient_pression")
952 {
953 const RefObjU& modele_turbulence = ch_inc_P1NC.equation().get_modele(TURBULENCE);
954 if (modele_turbulence && sub_type(Modele_turbulence_hyd_base,modele_turbulence.valeur()))
955 {
956 const Modele_turbulence_hyd_base& mod_turb = ref_cast(Modele_turbulence_hyd_base,modele_turbulence.valeur());
957 const Turbulence_paroi_base& loipar = mod_turb.loi_paroi();
958 DoubleTab tau_tan;
959
960 if (mod_turb.utiliser_loi_paroi())
961 {
962 tau_tan.ref(loipar.Cisaillement_paroi());
963
964 // Interpolation des vitesses a la paroi
965
966 const Domaine_VEF& domaine_VE_chF = ch_inc_P1NC.domaine_vef();
967 const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF,ch_inc_P1NC.equation().domaine_Cl_dis());
968 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
969 const Domaine& domchF = domaine_VE_chF.domaine();
970 const IntTab& face_voisins = domaine_VE_chF.face_voisins();
971 const IntTab& elem_faces = domaine_VE_chF.elem_faces();
972 int nfac = domchF.nb_faces_elem(0);
973 int nb_cl=les_cl.size();
974 IntVect num(nfac);
975 int num_cl,fac,ind_fac;
976
977 double kappa=0.41;
978 double norm_v,val0,val1,val2,norm_tau,u_tau;
979
980 for (num_cl=0; num_cl<nb_cl; num_cl++)
981 {
982 //Boucle sur les bords
983 const Cond_lim& la_cl = les_cl[num_cl];
984 const Front_VF& la_front_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
985 int ndeb = 0;
986 int nfin = la_front_dis.nb_faces_tot();
987
988 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) && tau_tan.size()!=0)
989 {
990 for (ind_fac=ndeb; ind_fac<nfin ; ind_fac++)
991 {
992 fac = la_front_dis.num_face(ind_fac);
993 int elem = face_voisins(fac,0);
994
995 if (Objet_U::dimension == 2 )
996 {
997 num[0]=elem_faces(elem,0);
998 num[1]=elem_faces(elem,1);
999 if (num[0]==fac) num[0]=elem_faces(elem,2);
1000 else if (num[1]==fac) num[1]=elem_faces(elem,2);
1001
1002 norm_v=norm_2D_vit1_lp(valeurs,fac,num[0],num[1],domaine_VE_chF,val0,val1);
1003 norm_tau=sqrt(tau_tan(fac,0)*tau_tan(fac,0)+tau_tan(fac,1)*tau_tan(fac,1));
1004 u_tau=sqrt(norm_tau);
1005 /*
1006 // PQ : 06/03 : pour un champ de vitesse 1D, l'incompressibilite genere des vitesses (meme faibles) sur toutes
1007 // les composantes, qui se repercutent par la suite au niveau de la CL ici.
1008 // (se manifestant en cours de calcul en entree par un pic de vitesse en proche paroi).
1009 // Pour assurer des entrees "correctes" (1D) on remedie a ce probleme via l'operation suivante
1010 // (en supposant que les entrees sont alignees avec une des directions d'espace)
1011
1012 if(std::fabs(val0)>0.999) {val0=1.;val1=0.;}
1013 if(std::fabs(val1)>0.999) {val0=0.;val1=1.;}
1014 */
1015
1016 valeurs(fac,0)=(norm_v-u_tau/kappa)*val0;
1017 valeurs(fac,1)=(norm_v-u_tau/kappa)*val1;
1018
1019 // for (i=0 ; i<Objet_U::dimension ; i++)
1020 // valeurs(fac,i)=(valeurs(num[0],i)+valeurs(num[1],i))/2.-sqrt(tau_tan(fac,i))/kappa;
1021 // valeurs(fac,0)=(valeurs(num[0],0)+valeurs(num[1],0))/2.-sqrt(tau_tan(fac,0))/kappa;
1022 }
1023 else
1024 {
1025 num[0]=elem_faces(elem,0);
1026 num[1]=elem_faces(elem,1);
1027 num[2]=elem_faces(elem,2);
1028 if (num[0]==fac) num[0]=elem_faces(elem,3);
1029 else if (num[1]==fac) num[1]=elem_faces(elem,3);
1030 else if (num[2]==fac) num[2]=elem_faces(elem,3);
1031
1032 norm_v=norm_3D_vit1_lp(valeurs,fac,num[0],num[1],num[2],domaine_VE_chF,val0,val1,val2);
1033 norm_tau=sqrt(tau_tan(fac,0)*tau_tan(fac,0)+tau_tan(fac,1)*tau_tan(fac,1)+tau_tan(fac,2)*tau_tan(fac,2));
1034 u_tau=sqrt(norm_tau);
1035 /*
1036 // PQ : 06/03 : pour un champ de vitesse 1D, l'incompressibilite genere des vitesses (meme faibles) sur toutes
1037 // les composantes, qui se repercutent par la suite au niveau de la CL ici.
1038 // (se manifestant en cours de calcul en entree par un pic de vitesse en proche paroi).
1039 // Pour assurer des entrees "correctes" (1D) on remedie a ce probleme via l'operation suivante
1040 // (en supposant que les entrees sont alignees suivant une des directions d'espace)
1041
1042 if(std::fabs(val0)>0.999) {val0=1.;val1=0.;val2=0.;}
1043 if(std::fabs(val1)>0.999) {val0=0.;val1=1.;val2=0.;}
1044 if(std::fabs(val2)>0.999) {val0=0.;val1=0.;val2=1.;}
1045 */
1046 valeurs(fac,0)=(norm_v-u_tau/kappa)*val0;
1047 valeurs(fac,1)=(norm_v-u_tau/kappa)*val1;
1048 valeurs(fac,2)=(norm_v-u_tau/kappa)*val2;
1049
1050
1051 // DoubleTab xv=domaine_VEF.xv();
1052 // Cerr << fac << " " << val0 << " " << val1 << " " << val2 <<finl;
1053 // Cerr << fac << " " << xv(fac,2) << " " << valeurs(fac,0) << " " << valeurs(fac,1) << " " << valeurs(fac,2) << " " << u_tau << " " << norm_v <<finl;
1054
1055 // for (i=0 ; i<Objet_U::dimension ; i++)
1056 // valeurs(fac,i)=(valeurs(num[0],i)+valeurs(num[1],i)+valeurs(num[2],i))/3.-sqrt(tau_tan(fac,i))/kappa;
1057 // valeurs(fac,0)=(valeurs(num[0],0)+valeurs(num[1],0)+valeurs(num[2],0))/3.-sqrt(tau_tan(fac,0))/kappa;
1058 }
1059 }//fac
1060 }//Dirichlet_paroi
1061 } //Boucle sur les bords
1062 valeurs.echange_espace_virtuel();
1063 }//!Paroi_negligeable_VEF
1064 }
1065 }
1066 }
1067
1068 DoubleTab valeurs_fac(cha.valeurs());
1069 cha.valeurs()=valeurs;
1070 Debog::verifier(" verif valeurs " , valeurs);
1071
1072 if (sub_type(Champ_P1NC,cha))
1073 {
1074 Champ_P1NC& ch_inc = ref_cast(Champ_P1NC,cha);
1075 valeur_P1_L2(ch_inc,dom);
1076 }
1077 else if (sub_type(Champ_Fonc_P1NC,cha))
1078 {
1079 Champ_Fonc_P1NC& ch_fonc = ref_cast(Champ_Fonc_P1NC,cha);
1080 valeur_P1_L2(ch_fonc,dom);
1081 }
1082
1083 Debog::verifier(" verif valeurs " , valeurs);
1084 cha.valeurs()=valeurs_fac;
1085 valeurs_fac=valeurs;
1086 Debog::verifier(" verif ch_som_ " , ch_som_);
1087 for(int face=0; face<nb_faces; face++)
1088 {
1089 int s1=domaine_VEF.face_sommets(face,0);
1090 s1=dom.get_renum_som_perio(s1);
1091 int s2=domaine_VEF.face_sommets(face,1);
1092 s2=dom.get_renum_som_perio(s2);
1093 int s3=-1;
1094 if(Objet_U::dimension==3)
1095 s3=dom.get_renum_som_perio(domaine_VEF.face_sommets(face,2));
1096 for(int comp=0; comp<nb_comp; comp++)
1097 {
1098 double somme=0;
1099 somme+=(ch_som_(s1,comp)+ch_som_(s2,comp));
1100 if(Objet_U::dimension==3)
1101 somme+=ch_som_(s3,comp);
1102 valeurs(face, comp)=coeff*somme;
1103 }
1104 }
1105 valeurs.echange_espace_virtuel();
1106 Debog::verifier(" verif valeurs L2" , valeurs);
1107}
1108
1109
1110void Champ_P1NC_implementation::filtrer_resu(DoubleTab& resu) const
1111{
1112 const Champ_P1NC& cha_const =ref_cast(Champ_P1NC,le_champ());
1113 Champ_P1NC& cha=ref_cast_non_const(Champ_P1NC, cha_const);
1114 double non_prepare=0;
1115 if (cha.MatP1NC2P1_L2.nb_lignes() <2)
1116 non_prepare=1;
1117 non_prepare=Process::mp_min(non_prepare);
1118 if (non_prepare==1)
1119 {
1120 // GF si on n a pas prepare la matrice, on appelle valeur_P1_L2(cha, cha.domaine()); et on remet les anciennes valeurs aux sommets
1121 DoubleTab vit_som_sa=cha.ch_som();
1122 valeur_P1_L2(cha, cha.domaine());
1123 cha.ch_som()=vit_som_sa;
1124 }
1125
1126 const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF,
1127 cha.equation().domaine_Cl_dis());
1128 const Domaine_VEF& domaine_VEF = domaine_vef();
1129 const DoubleVect& volumes_entrelaces= domaine_VEF.volumes_entrelaces();
1130 const Domaine& dom =domaine_VEF.domaine();
1131 const IntTab& face_sommets=domaine_VEF.face_sommets();
1132 int nb_faces=resu.dimension(0);
1133 int nb_faces_tot=resu.dimension_tot(0);
1134 int n_bord,nb_comp=resu.line_size();
1135
1136 DoubleTab resu_som(cha.ch_som());
1137 resu_som = 0;
1138
1139 ArrOfInt virtfait(nb_faces_tot-nb_faces);
1140
1141 double coeff= 1./Objet_U::dimension;
1142 int face,ind_face;
1143 {
1144 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
1145 {
1146 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1147 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1148 int nb_faces_tot_bord=le_bord.nb_faces_tot();
1149 int nb_faces_bord = le_bord.nb_faces();
1150 int num1 = 0;
1151 int num2 = nb_faces_tot_bord;
1152
1153 if (sub_type(Periodique,la_cl.valeur()))
1154 {
1155 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
1156 int face_associee;
1157 IntVect fait(nb_faces_tot_bord);
1158 for (ind_face=num1; ind_face<num2; ind_face++)
1159 {
1160 face = le_bord.num_face(ind_face);
1161 if (fait(ind_face) == 0)
1162 {
1163 fait(ind_face) = 1;
1164 if (ind_face-nb_faces_bord>=0) virtfait[face-nb_faces]=1;
1165 face_associee=la_cl_perio.face_associee(ind_face);
1166 fait(face_associee) = 1;
1167 if (ind_face-nb_faces_bord>=0)
1168 {
1169 face_associee = le_bord.num_face(face_associee);
1170 virtfait[face_associee-nb_faces]=1;
1171 }
1172
1173 for(int isom=0; isom<Objet_U::dimension; isom++)
1174 {
1175 int som=face_sommets(face,isom);
1176 som=dom.get_renum_som_perio(som);
1177 for (int comp=0; comp<nb_comp; comp++)
1178 resu_som(som, comp)+=coeff*resu(face, comp);
1179 }
1180
1181 }
1182 }
1183 }
1184 else
1185 {
1186 for (ind_face=num1; ind_face<num2; ind_face++)
1187 {
1188 face = le_bord.num_face(ind_face);
1189 if (ind_face-nb_faces_bord>=0) virtfait[face-nb_faces]=1;
1190 for(int isom=0; isom<Objet_U::dimension; isom++)
1191 {
1192 int som=face_sommets(face,isom);
1193 som=dom.get_renum_som_perio(som);
1194 for (int comp=0; comp<nb_comp; comp++)
1195 resu_som(som, comp)+=coeff*resu(face, comp);
1196 }
1197 }
1198 }
1199 }
1200
1201
1202 int deb=domaine_VEF.premiere_face_int();
1203
1204
1205 for(face=deb; face<nb_faces_tot; face++)
1206 {
1207 int indvirtface=face-nb_faces;
1208 // on regarde si la face virtuelle n'a pas deja ete traitee
1209 if (indvirtface>=0)
1210 if (virtfait[indvirtface]==1) continue;
1211 for(int isom=0; isom<Objet_U::dimension; isom++)
1212 {
1213 int som=face_sommets(face,isom);
1214 som=dom.get_renum_som_perio(som);
1215 for(int comp=0; comp<nb_comp; comp++)
1216 resu_som(som, comp)+=coeff*resu(face, comp);
1217 }
1218 }
1219 }
1220
1221 // resu_som.echange_espace_virtuel();
1222 // Debog::verifier("Champ_P1NC_implementation::filtrer_resu 1.5 ",resu);
1223 // Debog::verifier("Champ_P1NC_implementation::filtrer_resu_som 1.5 ",resu_som);
1224
1225 // ************************************************************************************************
1226 // ************************************************************************************************
1227 //
1228 // Pour imposer residu_som = residu donne par CL
1229 //
1230 // ************************************************************************************************
1231 // ************************************************************************************************
1232 IntVect test_cl_imposee(resu_som.dimension_tot(0));
1233 test_cl_imposee = 0;
1234
1235 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
1236 {
1237 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1238 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1239 int num1 = 0;
1240 int num2 = le_bord.nb_faces_tot() ;
1241
1242 // if (sub_type(Dirichlet,la_cl.valeur()) || sub_type(Dirichlet_homogene,la_cl.valeur()))
1243 if (sub_type(Dirichlet_homogene,la_cl.valeur()))
1244 {
1245 for (ind_face=num1; ind_face<num2; ind_face++)
1246 {
1247 face = le_bord.num_face(ind_face);
1248 double vol=volumes_entrelaces(face);
1249 for(int isom=0; isom<Objet_U::dimension; isom++)
1250 {
1251 int som=face_sommets(face,isom);
1252 som=dom.get_renum_som_perio(som);
1253 for (int comp=0; comp<nb_comp; comp++)
1254 {
1255 if (test_cl_imposee(som)==0)
1256 resu_som(som,comp) = 0;
1257 resu_som(som, comp) += coeff_penalisation * resu(face, comp)/vol;
1258 }
1259 test_cl_imposee(som) += 1;
1260 }
1261 }
1262 }
1263 }
1264
1265 resu_som.echange_espace_virtuel();
1266
1267 SolveurSys& solv=ref_cast(SolveurSys,cha.solveur_L2);
1268 int i,nbs = resu_som.dimension(0);
1269
1270 DoubleVect secmem(cha.ch_som_vect());
1271 DoubleVect solution(cha.ch_som_vect());
1272
1273 // Modif CD 15/04/03 pour paralleliser le filtrage
1274 // Matrice matrice_tmp;
1275 // cha.dimensionner_Mat_Bloc_Morse_Sym(matrice_tmp);
1276 // cha.Mat_Morse_to_Mat_Bloc(matrice_tmp);
1277 // Fin Modif CD 15/04/03
1278
1279 // Debog::verifier_Mat_sommet("Champ_P1NC_implementation::filtrer_resu cha.MatP1NC2P1_L2 ", cha.MatP1NC2P1_L2);
1280
1281 for(int k=0; k<nb_comp; k++)
1282 {
1283 solution=0.;
1284 for(i=0; i<nbs; i++)
1285 secmem(i)=resu_som(i,k);
1286
1287 secmem.echange_espace_virtuel();
1288
1289 solv.resoudre_systeme(cha.MatP1NC2P1_L2_Parallele.valeur(), secmem, solution);
1290
1291 for(i=0; i<nbs; i++)
1292 resu_som(i,k)=solution(i);
1293 }
1294
1295 resu=0;
1296 {
1297 for(face=0; face<nb_faces; face++)
1298 {
1299 for(int isom=0; isom<Objet_U::dimension; isom++)
1300 {
1301 int som=dom.get_renum_som_perio(domaine_VEF.face_sommets(face,isom));
1302 for(int comp=0; comp<nb_comp; comp++)
1303 resu(face, comp)+=coeff*resu_som(som, comp);
1304 }
1305 }
1306 }
1307 {
1308 for(face=0; face<nb_faces; face++)
1309 {
1310 double vol=volumes_entrelaces(face);
1311 for(int comp=0; comp<nb_comp; comp++)
1312 resu(face, comp)*=vol;
1313 }
1314 }
1315}
1316
1317DoubleVect& Champ_P1NC_implementation::valeur_a_elem(const DoubleVect& position,
1318 DoubleVect& val,
1319 int le_poly) const
1320{
1321 const Champ_base& cha=le_champ();
1322 const int N = cha.nb_comp(), D = Objet_U::dimension;
1323 const Domaine& domaine_geom = get_domaine_geom();
1324 const Domaine_VEF& domaine_VEF = domaine_vef();
1325 const DoubleTab& coord = domaine_geom.coord_sommets();
1326 const IntTab& sommet_poly = domaine_geom.les_elems();
1327 const IntTab& elem_faces=domaine_VEF.elem_faces();
1328 const DoubleTab& ch = cha.valeurs();
1329 val = 0;
1330
1331 if (le_poly != -1)
1332 {
1333 const double xs = position(0), ys = position(1),
1334 zs = (D == 3) ? position(2) : 0;
1335 for (int i = 0; i < D + 1; i++)
1336 {
1337 const int f = elem_faces(le_poly, i);
1338 for(int n = 0; n < N; n++)
1339 val(n) += ch(f, n) * ((D == 2) ? fonction_forme_2D(xs, ys, le_poly, i, sommet_poly, coord) : fonction_forme_3D(xs, ys, zs, le_poly, i, sommet_poly, coord));
1340 }
1341 }
1342 return val;
1343}
1344
1345double Champ_P1NC_implementation::
1346valeur_a_elem_compo(const DoubleVect& position, int le_poly, int ncomp) const
1347{
1348 //Cerr << "Champ_P1NC_implementation::valeur_a_elem_compo" << finl;
1349 const Champ_base& cha=le_champ();
1350 const Domaine& domaine_geom = get_domaine_geom();
1351 const int D = Objet_U::dimension;
1352 const Domaine_VEF& domaine_VEF = domaine_vef();
1353 const DoubleTab& coord = domaine_geom.coord_sommets();
1354 const IntTab& sommet_poly = domaine_geom.les_elems();
1355 const IntTab& elem_faces=domaine_VEF.elem_faces();
1356 const DoubleTab& ch = cha.valeurs();
1357 double val = 0.;
1358
1359 if (le_poly != -1)
1360 {
1361 // Calcul d'apres les fonctions de forme sur le triangle ou le tetraedre
1362 const double xs = position(0), ys = position(1),
1363 zs = (D == 3) ? position(2) : 0;
1364 for (int i = 0; i < D + 1; i++)
1365 {
1366 const int f = elem_faces(le_poly,i);
1367 val += ch(f, ncomp) * ((D == 2) ? fonction_forme_2D(xs, ys, le_poly, i, sommet_poly, coord) : fonction_forme_3D(xs, ys, zs, le_poly, i, sommet_poly, coord));
1368 }
1369 }
1370
1371 return val;
1372}
1373
1374/**
1375* @brief Computes values at the centers of gravity for a P1NC field
1376 *
1377 * This method calculates the values of a P1NC field at the centers of gravity
1378 * of the mesh elements. For each component of the field, it averages the values
1379 * from the faces of each polygon.
1380*
1381* @param[in,out] val DoubleTab that will store the computed values.
1382* Must have 1 or 2 dimensions.
1383* If 1D, it will be resized to (dimension_tot(0), 1).
1384*
1385* @return Reference to the modified val DoubleTab containing computed values
1386*
1387* @throws Process::exit() if val has more than 2 dimensions
1388*
1389* @note Implementation uses GPU parallelization through Kokkos
1390* @note For each polygon, the value is computed as average of D+1 face values,
1391 * where D is the spatial dimension
1392*/
1393DoubleTab& Champ_P1NC_implementation::valeur_aux_centres_de_gravite(const Domaine& dom, DoubleTab& tab_val) const
1394{
1395 const Champ_base& cha = le_champ();
1396 int nb_compo_ = cha.nb_comp();
1397 const Domaine_VEF& domaine_VEF = domaine_vef();
1398 if (domaine_VEF.domaine() != dom) Process::exit("Error, you must use valeur_aux_centres_de_gravite() on the whole discretized mesh.");
1399 if (tab_val.nb_dim() > 2)
1400 {
1401 Cerr << "Erreur TRUST dans Champ_P1NC_implementation::valeur_aux_elems()\n";
1402 Cerr << "Le DoubleTab val a plus de 2 entrees\n";
1403 Process::exit();
1404 }
1405
1406 // TODO : FIXME
1407 // For FT the resize should be done in its good position and not here ...
1408 if (tab_val.nb_dim() == 1) tab_val.resize(tab_val.dimension_tot(0), 1);
1409
1410 int D = Objet_U::dimension;
1411 int nb_elem = tab_val.dimension(0);
1412 CDoubleTabView ch = cha.valeurs().view_ro();
1413 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
1414 DoubleTabView val = tab_val.view_wo();
1415 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0,0}, {nb_elem,nb_compo_}), KOKKOS_LAMBDA(const int le_poly, const int ncomp)
1416 {
1417 double sum = 0;
1418 for (int i = 0; i < D + 1; i++)
1419 {
1420 int face = elem_faces(le_poly, i);
1421 sum += ch(face, ncomp);
1422 }
1423 val(le_poly, ncomp) = sum / (D + 1);
1424 });
1425 end_gpu_timer(__KERNEL_NAME__);
1426 return tab_val;
1427}
1428
1429DoubleTab& Champ_P1NC_implementation::
1430valeur_aux_elems(const DoubleTab& positions,
1431 const IntVect& les_polys,
1432 DoubleTab& tab_val) const
1433{
1434 if (les_polys.size() == 0) return tab_val;
1435 const Domaine_VEF& domaine_VEF = domaine_vef();
1436#ifdef TRUST_USE_GPU
1437 if (les_polys.size()==domaine_VEF.nb_elem())
1438 Cerr << "Warning, Champ_P1NC_implementation::valeur_aux_elems(...) called whereas call to Champ_P1NC_implementation::valeur_aux_centres_de_gravite(...) could be more efficient!" << finl;
1439#endif
1440 const Champ_base& cha=le_champ();
1441 int nb_compo_=cha.nb_comp();
1442 const DoubleTab& ch = cha.valeurs();
1443
1444 const Domaine& domaine_geom = get_domaine_geom();
1445 const DoubleTab& coord = domaine_geom.coord_sommets();
1446 const IntTab& sommet_poly = domaine_geom.les_elems();
1447 if (tab_val.nb_dim() > 2)
1448 {
1449 Cerr << "Erreur TRUST dans Champ_P1NC_implementation::valeur_aux_elems()\n";
1450 Cerr << "Le DoubleTab val a plus de 2 entrees\n";
1451 Process::exit();
1452 }
1453
1454 // TODO : FIXME
1455 // For FT the resize should be done in its good position and not here ...
1456 if (tab_val.nb_dim() == 1) tab_val.resize(tab_val.dimension_tot(0),1);
1457
1458 int D = Objet_U::dimension;
1459 CIntTabView sommet_poly_v = sommet_poly.view_ro();
1460 CDoubleTabView coord_v = coord.view_ro();
1461 CDoubleTabView ch_v = ch.view_ro();
1462 CIntTabView elem_faces_v = domaine_VEF.elem_faces().view_ro();
1463 CDoubleTabView positions_v = positions.view_ro();
1464 CIntArrView les_polys_v = les_polys.view_ro();
1465 DoubleTabView val = tab_val.view_wo();
1466 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0,0}, {les_polys.size(),nb_compo_}), KOKKOS_LAMBDA(
1467 const int rang_poly, const int ncomp)
1468 {
1469 int le_poly=les_polys_v(rang_poly);
1470 if (le_poly == -1)
1471 val(rang_poly, ncomp) = 0;
1472 else
1473 {
1474 double xs = positions_v(rang_poly,0);
1475 double ys = positions_v(rang_poly,1);
1476 double zs = (D == 3) ? positions_v(rang_poly,2) : 0;
1477 double sum = 0;
1478 for (int i = 0; i < D + 1; i++)
1479 {
1480 int face = elem_faces_v(le_poly, i);
1481 sum += ch_v(face, ncomp) *
1482 ((D == 2) ? fonction_forme_2D_v(xs, ys, le_poly, i, sommet_poly_v, coord_v) : fonction_forme_3D_v(xs, ys, zs, le_poly, i, sommet_poly_v, coord_v));
1483 }
1484 val(rang_poly, ncomp) = sum;
1485 }
1486 });
1487 end_gpu_timer(__KERNEL_NAME__);
1488 return tab_val;
1489}
1490
1491DoubleVect& Champ_P1NC_implementation::
1492valeur_aux_elems_compo(const DoubleTab& positions,
1493 const IntVect& les_polys,
1494 DoubleVect& val,int ncomp) const
1495{
1496 if (les_polys.size() == 0) return val;
1497 int les_polys_size = les_polys.size();
1498 const Champ_base& cha=le_champ();
1499 double xs,ys,zs;
1500 int face;
1501 const Domaine_VEF& domaine_VEF = domaine_vef();
1502 const Domaine& domaine_geom = get_domaine_geom();
1503 const DoubleTab& coord = domaine_geom.coord_sommets();
1504 const IntTab& sommet_poly = domaine_geom.les_elems();
1505 // Commenter en attendant de comprendre le dimensionnement de xp apres la sortie de version
1506 // assert(val.size() == les_polys_size);
1507 int le_poly, D = Objet_U::dimension;
1508
1509 const DoubleTab& ch = cha.valeurs();
1510 ToDo_Kokkos("critical");
1511 for(int rang_poly=0; rang_poly<les_polys_size; rang_poly++)
1512 {
1513 le_poly=les_polys(rang_poly);
1514 if (le_poly == -1) val(rang_poly) = 0;
1515 else
1516 {
1517 // Calcul d'apres les fonctions de forme sur le triangle ou le tetraedre
1518 val(rang_poly) = 0;
1519 xs = positions(rang_poly,0);
1520 ys = positions(rang_poly,1);
1521 zs = (D == 3) ? positions(rang_poly,2) : 0.;
1522 for (int i=0; i < D + 1; i++)
1523 {
1524 face = domaine_VEF.elem_faces(le_poly,i);
1525 val(rang_poly) += ch(face, ncomp) *
1526 ( (D == 2) ? fonction_forme_2D(xs,ys,le_poly,i,sommet_poly,coord) :
1527 fonction_forme_3D(xs,ys,zs,le_poly,i,sommet_poly,coord));
1528 }
1529 }
1530 }
1531
1532 return val;
1533}
1534
1535
1536DoubleTab& Champ_P1NC_implementation::
1537valeur_aux_elems_smooth(const DoubleTab& positions,
1538 const IntVect& les_polys,
1539 DoubleTab& val)
1540{
1541 if ((!sub_type(Champ_P1NC,le_champ()))&&(!sub_type(Champ_Fonc_P1NC,le_champ())))
1542 {
1543 Cerr << "L'option chsom des sondes ne s'applique pour le moment que pour les Champ_P1NC en VEF" << finl;
1544 Cerr << "or votre champ " << le_champ().le_nom() << " est de type " << le_champ().que_suis_je() << "." << finl;
1545 // Je rajoute le message ci-dessous car le message habituel n'apparait pas !
1546 Cerr << "TRUST a provoque une erreur et va s'arreter" << finl;
1547 Process::exit();
1548 }
1549
1550 Champ_base& cha =ref_cast(Champ_base,le_champ());
1551 const int les_polys_size = les_polys.size(), nb_compo_=cha.nb_comp(), D = Objet_U::dimension;
1552 const DoubleTab& ch_sommet= (ch_som());
1553 const Domaine_VEF& domaine_VEF = domaine_vef();
1554 const Domaine& domaine_geom = get_domaine_geom();
1555 const Domaine& dom=domaine_VEF.domaine();
1556 const DoubleTab& coord = domaine_geom.coord_sommets();
1557 const IntTab& sommet_poly = domaine_geom.les_elems();
1558 if (val.nb_dim() == 2)
1559 {
1560 assert((val.dimension(0) == les_polys_size)||(val.dimension_tot(0) == les_polys_size));
1561 assert(val.line_size() == nb_compo_);
1562 }
1563 else
1564 {
1565 Cerr << "Erreur TRUST dans Champ_P1NC_implementation::valeur_aux_elems_smooth()\n";
1566 Cerr << "Le DoubleTab val a plus de 2 entrees\n";
1567 Process::exit();
1568 }
1569 if (!filtrer_L2_deja_appele_)
1570 {
1571 // Filtrer L2 ne marche que pour les vecteurs
1572 // C.MALOD 19/12/2006 : Ce n'est plus vrai, ca marche aussi avec les scalaires.
1573 if (nb_compo_>0)
1574 {
1575 DoubleTab val_sauv=cha.valeurs();
1576 filtrer_L2(cha.valeurs());
1577 cha.valeurs()=val_sauv;
1578 }
1579 filtrer_L2_deja_appele_=1;
1580 }
1581
1582 // calcul de la valeur aux elements suivant les coordonnees barycentriques (repris de Champ_P1)
1583 val = 0.;
1584 int p;
1585 ToDo_Kokkos("critical");
1586 for(int rang_poly=0; rang_poly<les_polys_size; rang_poly++)
1587 if ((p = les_polys(rang_poly)) != -1)
1588 {
1589 const double xs = positions(rang_poly,0), ys = positions(rang_poly,1),
1590 zs = (D == 3) ? positions(rang_poly,2) : 0.;
1591 for(int ncomp=0; ncomp<nb_compo_; ncomp++)
1592 for (int i = 0; i < D + 1; i++)
1593 {
1594 int som = sommet_poly(p, i);
1595 int som1=dom.get_renum_som_perio(som);
1596 val(rang_poly, ncomp) += ch_sommet(som1, ncomp) * ((D == 2) ? coord_barycentrique(sommet_poly, coord, xs, ys, p, i) : coord_barycentrique(sommet_poly, coord, xs, ys, zs, p, i));
1597 }
1598 }
1599 return val;
1600}
1601
1602
1603DoubleVect& Champ_P1NC_implementation::
1604valeur_aux_elems_compo_smooth(const DoubleTab& positions,
1605 const IntVect& les_polys,
1606 DoubleVect& val,int ncomp)
1607{
1608 if ((!sub_type(Champ_P1NC,le_champ()))&&(!sub_type(Champ_Fonc_P1NC,le_champ())))
1609 {
1610 Cerr << "L'option chsom des sondes ne s'applique pour le moment que pour les Champ_P1NC en VEF" << finl;
1611 Cerr << "or votre champ " << le_champ().le_nom() << " est de type " << le_champ().que_suis_je() << "." << finl;
1612 // Je rajoute le message ci-dessous car le message habituel n'apparait pas !
1613 Cerr << "TRUST a provoque une erreur et va s'arreter" << finl;
1614 Process::exit();
1615 }
1616 Champ_base& cha =le_champ();
1617 const int les_polys_size = les_polys.size(), nb_compo_=cha.nb_comp(), D = Objet_U::dimension;
1618 const DoubleTab& ch_sommet= (ch_som());
1619 const Domaine_VEF& domaine_VEF = domaine_vef();
1620 const Domaine& dom=domaine_VEF.domaine();
1621 const Domaine& domaine_geom = get_domaine_geom();
1622 const DoubleTab& coord = domaine_geom.coord_sommets();
1623 const IntTab& sommet_poly = domaine_geom.les_elems();
1624 if (!filtrer_L2_deja_appele_)
1625 {
1626 // Filtrer L2 ne marche que pour les vecteurs
1627 // C.MALOD 19/12/2006 : Ce n'est plus vrai, ca marche aussi avec les scalaires.
1628 if (nb_compo_>0)
1629 {
1630 DoubleTab val_sauv=cha.valeurs();
1631 filtrer_L2(cha.valeurs());
1632 cha.valeurs()=val_sauv;
1633 }
1634 filtrer_L2_deja_appele_=1;
1635 }
1636 assert(val.size_totale() >= les_polys_size);
1637 val = 0.;
1638 int p;
1639
1640 // calcul de la valeur aux elements suivant les coordonnees barycentriques (repris de Champ_P1)
1641 ToDo_Kokkos("critical");
1642 for(int rang_poly=0; rang_poly<les_polys_size; rang_poly++)
1643 if ((p = les_polys(rang_poly)) != -1)
1644 {
1645 const double xs = positions(rang_poly,0), ys = positions(rang_poly,1),
1646 zs = (D == 3) ? positions(rang_poly,2) : 0.;
1647 for (int i = 0; i < D + 1; i++)
1648 {
1649 int som = sommet_poly(p, i);
1650 int som1=dom.get_renum_som_perio(som);
1651 val(rang_poly) += ch_sommet(som1, ncomp)
1652 * ((D == 2) ? coord_barycentrique(sommet_poly, coord, xs, ys, p, i) : coord_barycentrique(sommet_poly, coord, xs, ys, zs, p, i));
1653 }
1654 }
1655 return val;
1656}
1657
1658DoubleTab& Champ_P1NC_implementation::
1659valeur_aux_sommets(const Domaine& dom,
1660 DoubleTab& tab_champ_som) const
1661{
1662 int nb_som = dom.nb_som_tot();
1663 int nb_som_return = tab_champ_som.dimension_tot(0);
1664 int nb_compo_ = le_champ().nb_comp();
1665 IntTrav tab_compteur(nb_som);
1666
1667 const Domaine_VEF& domaine_VEF = domaine_vef();
1668 assert(domaine_vef().domaine()==dom);
1669 if (methode_calcul_valeurs_sommets_static>0)
1670 {
1671 int nb_elem_tot = dom.nb_elem_tot();
1672 int nb_som_elem = dom.nb_som_elem();
1673
1674 DoubleTrav tab_min_som(nb_som,nb_compo_) ;
1675 DoubleTrav tab_max_som(nb_som,nb_compo_) ;
1676
1677 tab_min_som = 1.e+30;
1678 tab_max_som = 1.e-30; // ToDo bug ? tab_max_som = -1.e+30
1679 if (getenv("TRUST_FIX_P1NC_POST_SOM")!=nullptr) tab_max_som = -1.e+30;
1680 tab_compteur = 0;
1681 tab_champ_som = 0.;
1682 // Loop on elem (~30ms -> ~8ms by using parallelism on the 2 nested loops with MDRangePolicy)
1683 {
1684 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
1685 CIntTabView sommet_elem = dom.les_elems().view_ro();
1686 CDoubleTabView ch = le_champ().valeurs().view_ro();
1687 IntArrView compteur = static_cast<ArrOfInt&>(tab_compteur).view_rw();
1688 DoubleTabView champ_som = tab_champ_som.view_rw();
1689 DoubleTabView min_som = tab_min_som.view_rw();
1690 DoubleTabView max_som = tab_max_som.view_rw();
1691 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0,0}, {nb_elem_tot,nb_som_elem}), KOKKOS_LAMBDA(
1692 const int num_elem, const int j)
1693 {
1694 int num_som = sommet_elem(num_elem, j);
1695 if (num_som < nb_som_return)
1696 {
1697 Kokkos::atomic_inc(&compteur(num_som));
1698 for (int ncomp = 0; ncomp < nb_compo_; ncomp++)
1699 {
1700 Kokkos::atomic_add(&champ_som(num_som, ncomp), valeur_a_sommet_compo(num_som, num_elem, ncomp, elem_faces, sommet_elem, ch));
1701 for (int face_adj = 0; face_adj < nb_som_elem; face_adj++)
1702 {
1703 if (face_adj != j)
1704 {
1705 int face_glob = elem_faces(num_elem, face_adj);
1706 Kokkos::atomic_min(&min_som(num_som, ncomp), ch(face_glob, ncomp));
1707 Kokkos::atomic_max(&max_som(num_som, ncomp), ch(face_glob, ncomp));
1708 }
1709 }
1710 }
1711 }
1712 });
1713 end_gpu_timer(__KERNEL_NAME__);
1714 }
1715 // Loop on nodes (~0.2ms, no gain by collapsing the 2 nested loops...)
1716 {
1717 int methode_calcul_valeurs_sommets = methode_calcul_valeurs_sommets_static; // Kokkos kernels can't use static var on device
1718 CDoubleTabView min_som = tab_min_som.view_ro();
1719 CDoubleTabView max_som = tab_max_som.view_ro();
1720 CIntArrView compteur = static_cast<ArrOfInt&>(tab_compteur).view_ro();
1721 DoubleTabView champ_som = tab_champ_som.view_rw();
1722 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0,0}, {nb_som_return,nb_compo_}), KOKKOS_LAMBDA(
1723 const int num_som, const int ncomp)
1724 {
1725 if (compteur(num_som)>0)
1726 {
1727 champ_som(num_som, ncomp) /= compteur(num_som);
1728 if (methode_calcul_valeurs_sommets > 1)
1729 {
1730 if (champ_som(num_som, ncomp) < min_som(num_som, ncomp))
1731 champ_som(num_som, ncomp) = min_som(num_som, ncomp);
1732 if (champ_som(num_som, ncomp) > max_som(num_som, ncomp))
1733 champ_som(num_som, ncomp) = max_som(num_som, ncomp);
1734 }
1735 }
1736 });
1737 end_gpu_timer(__KERNEL_NAME__);
1738 }
1739 }
1740 else
1741 {
1742 assert(methode_calcul_valeurs_sommets_static==-1);
1743 tab_champ_som = 0;
1744 tab_compteur = 0;
1745 int nb_comp=nb_compo_;
1746 int nb_faces_tot=domaine_VEF.nb_faces_tot();
1747 int nb_som_face=domaine_VEF.nb_som_face();
1748 const IntTab& face_sommets=domaine_VEF.face_sommets();
1749 const DoubleTab& ch = le_champ().valeurs();
1750 ToDo_Kokkos("critical"); // but never used !
1751 for(int face=0; face<nb_faces_tot; face++)
1752 {
1753
1754 for(int isom=0; isom<nb_som_face; isom++)
1755 {
1756 int som1=face_sommets(face,isom);
1757 if (som1<nb_som_return)
1758 {
1759 tab_compteur[som1]++;
1760 for(int k=0; k<nb_comp; k++)
1761 tab_champ_som(som1,k)+=ch(face,k);
1762 }
1763 }
1764 }
1765
1766
1767 for (int num_som=0; num_som<nb_som_return; num_som++)
1768 for (int ncomp=0; ncomp<nb_compo_; ncomp++)
1769 {
1770 tab_champ_som(num_som,ncomp) /= tab_compteur[num_som];
1771 }
1772 }
1773 return tab_champ_som;
1774}
1775
1776double Champ_P1NC_implementation::valeur_a_sommet_compo(int num_som, int num_elem, int ncomp) const
1777{
1778 const IntTab& elem_faces = domaine_vef().elem_faces();
1779 const IntTab& sommet_elem = get_domaine_geom().les_elems();
1780 const DoubleTab& ch = le_champ().valeurs();
1781 double val=0;
1782 if (num_elem != -1)
1783 {
1784 for (int i=0; i< Objet_U::dimension+1; i++)
1785 {
1786 int face = elem_faces(num_elem,i);
1787 if (sommet_elem(num_elem,i)==num_som)
1788 val -= (Objet_U::dimension-1)*ch(face, ncomp);
1789 else
1790 val += ch(face, ncomp);
1791 }
1792 }
1793 return val;
1794}
1795
1796KOKKOS_INLINE_FUNCTION
1797double Champ_P1NC_implementation::valeur_a_sommet_compo(int num_som, int num_elem, int ncomp, CIntTabView elem_faces, CIntTabView sommet_elem, CDoubleTabView ch) const
1798{
1799 double val=0;
1800 if (num_elem != -1)
1801 {
1802 int nb_som = (int)sommet_elem.extent(1);
1803 for (int i=0; i<nb_som; i++)
1804 {
1805 int face = elem_faces(num_elem,i);
1806 if (sommet_elem(num_elem,i)==num_som)
1807 val -= (nb_som-2)*ch(face, ncomp);
1808 else
1809 val += ch(face, ncomp);
1810 }
1811 }
1812 return val;
1813}
1814
1815//DoubleTab& Champ_P1NC_implementation::
1816//valeur_aux_sommets(const Domaine& dom,
1817// DoubleTab& ch_som) const
1818//{
1819// //Cerr << "Champ_P1NC_implementation::valeur_aux_sommets" << finl;
1820// Champ_P1NC& ch=ref_cast(Champ_P1NC, le_champ());
1821// valeur_P1_L2(ch, dom);
1822// DoubleTab& champ_sommet=ch.ch_som();
1823// if (champ_sommet.nb_dim()==1)
1824// {
1825// for (int k=0;k<champ_sommet.dimension(0);k++)
1826// ch_som(k,0)=champ_sommet(k);
1827// }
1828// else
1829// {
1830// for (int k=0;k<champ_sommet.dimension(0);k++)
1831// for (int ncomp=0;ncomp<champ_sommet.dimension(1);ncomp++)
1832// ch_som(k,ncomp)=champ_sommet(k,ncomp);
1833// }
1834// return ch_som;
1835//}
1836
1837DoubleVect& Champ_P1NC_implementation::
1838valeur_aux_sommets_compo(const Domaine& dom,
1839 DoubleVect& champ_som,
1840 int ncomp) const
1841{
1842 const DoubleTab& ch = le_champ().valeurs();
1843
1844 const Domaine_VEF& domaine_VEF = domaine_vef();
1845
1846
1847 int nb_som = dom.nb_som();
1848
1849 IntVect compteur(nb_som);
1850 if (methode_calcul_valeurs_sommets_static>0)
1851 {
1852 int nb_elem_tot = dom.nb_elem_tot();
1853 int nb_som_elem = dom.nb_som_elem();
1854 int num_elem,num_som,j;
1855 champ_som = 0;
1856 compteur = 0;
1857 const IntTab& elem_faces = domaine_VEF.elem_faces();
1858
1859
1860 DoubleVect min_som(nb_som) ;
1861 DoubleVect max_som(nb_som) ;
1862 int face_adj,face_glob;
1863
1864 for (j=0; j<nb_som; j++)
1865 {
1866 min_som(j) = 1.e+30 ;
1867 max_som(j) = 1.e-30 ;
1868 }
1869 ToDo_Kokkos("critical");
1870 for (num_elem=0; num_elem<nb_elem_tot; num_elem++)
1871 {
1872 for (j=0; j<nb_som_elem; j++)
1873 {
1874 num_som = dom.sommet_elem(num_elem,j);
1875 if(num_som < nb_som)
1876 {
1877 compteur[num_som]++;
1878
1879 champ_som(num_som) += valeur_a_sommet_compo(num_som,num_elem,ncomp);
1880
1881 for(face_adj=0; face_adj<nb_som_elem; face_adj ++)
1882 {
1883 if (face_adj != j )
1884 {
1885 face_glob = elem_faces(num_elem, face_adj);
1886 min_som(num_som) = std::min
1887 ( ch(face_glob,ncomp),min_som(num_som) ) ;
1888 max_som(num_som) = std::max
1889 ( ch(face_glob,ncomp),max_som(num_som) ) ;
1890 }
1891 }
1892
1893 }
1894 }
1895 }
1896 ToDo_Kokkos("critical");
1897 for (num_som=0; num_som<nb_som; num_som++)
1898 {
1899 champ_som(num_som) /= compteur[num_som];
1900 if (methode_calcul_valeurs_sommets_static>1)
1901 {
1902 if ( champ_som(num_som) < min_som(num_som) ) champ_som(num_som) = min_som(num_som) ;
1903 if ( champ_som(num_som) > max_som(num_som) ) champ_som(num_som) = max_som(num_som) ;
1904 }
1905
1906 }
1907 }
1908 else
1909 {
1910 champ_som = 0;
1911 compteur = 0;
1912 int nb_faces_tot=domaine_VEF.nb_faces_tot();
1913 int nb_som_face=domaine_VEF.nb_som_face();
1914 const IntTab& face_sommets=domaine_VEF.face_sommets();
1915 int face;
1916 ToDo_Kokkos("critical");
1917 for(face=0; face<nb_faces_tot; face++)
1918 {
1919
1920 for(int isom=0; isom<nb_som_face; isom++)
1921 {
1922 int som1=face_sommets(face,isom);
1923 if (som1<nb_som)
1924 {
1925 compteur[som1]++;
1926 champ_som(som1)+=ch(face,ncomp);
1927 }
1928 }
1929 }
1930
1931 ToDo_Kokkos("critical");
1932 for (int num_som=0; num_som<nb_som; num_som++)
1933 {
1934 champ_som(num_som) /= compteur[num_som];
1935 }
1936
1937 }
1938 return champ_som;
1939}
1940
1941DoubleTab& Champ_P1NC_implementation::remplir_coord_noeuds(DoubleTab& noeuds) const
1942{
1943 const Domaine_VEF& zvef = domaine_vef();
1944 const DoubleTab& xv = zvef.xv();
1945 int nb_fac = zvef.nb_faces();
1946 if ( (xv.dimension(0) == nb_fac ) && (xv.dimension(1) == Objet_U::dimension) )
1947 noeuds.ref(xv);
1948 else
1949 {
1950 Cerr << "Erreur dans Champ_P1NC_implementation::remplir_coord_noeuds()" << finl;
1951 Cerr << "Les centres de gravite des faces n'ont pas ete calcules" << finl;
1952 Process::exit();
1953 }
1954 return noeuds;
1955}
1956
1957
1958int Champ_P1NC_implementation::remplir_coord_noeuds_et_polys(DoubleTab& positions,
1959 IntVect& polys) const
1960{
1961 remplir_coord_noeuds(positions);
1962 const Domaine_VEF& zvef=domaine_vef();
1963 const IntTab& face_voisins=zvef.face_voisins();
1964 int nb_faces=zvef.nb_faces();
1965 polys.resize(nb_faces);
1966 ToDo_Kokkos("critical");
1967 for(int i= 0; i< nb_faces; i++)
1968 {
1969 polys(i)=face_voisins(i, 0);
1970 }
1971 return 0;
1972}
1973
1974int Champ_P1NC_implementation::imprime_P1NC(Sortie& os, int ncomp) const
1975{
1976 const Domaine_VEF& domaine_VEF=domaine_vef();
1977 const DoubleTab& pos_face=domaine_VEF.xv();
1978 const int nb_faces = domaine_VEF.nb_faces();
1979 const Champ_base& cha=le_champ();
1980 const DoubleTab& val = cha.valeurs();
1981 int fac;
1982 os << nb_faces << finl;
1983 for (fac=0; fac<nb_faces; fac++)
1984 {
1985 if (Objet_U::dimension==3)
1986 os << pos_face(fac,0) << " " << pos_face(fac,1) << " " << pos_face(fac,2) << " " ;
1987 if (Objet_U::dimension==2)
1988 os << pos_face(fac,0) << " " << pos_face(fac,1) << " " ;
1989 os << val(fac,ncomp) << finl;
1990 }
1991 os << finl;
1992 Cout << "Champ_P1NC_implementation::imprime_P1NC FIN >>>>>>>>>> " << finl;
1993 return 1;
1994}
1995
1996int Champ_P1NC_implementation::fixer_nb_valeurs_nodales(int n)
1997{
1998 ch_som_.reset();
1999 const int nb_dim = le_champ().valeurs().nb_dim();
2000 if (nb_dim > 1)
2001 ch_som_.resize(0, le_champ().nb_comp());
2002 const Domaine& dom = domaine_vef().domaine();
2003 dom.creer_tableau_sommets(ch_som_);
2004 ch_som_vect_.reset();
2005 dom.creer_tableau_sommets(ch_som_vect_);
2006 return n;
2007}
2008// Ajout de methodes pour le filtrage L2 en //
2009int Champ_P1NC_implementation::nb_colonnes_tot()
2010{
2011 if (sub_type(Champ_P1NC,le_champ()))
2012 {
2013 Champ_P1NC& cha =ref_cast(Champ_P1NC,le_champ());
2014
2015 return cha.ch_som_vect().size_totale();
2016 }
2017 else if (sub_type(Champ_Fonc_P1NC,le_champ()))
2018 {
2019 Champ_Fonc_P1NC& cha = ref_cast(Champ_Fonc_P1NC,le_champ());
2020 return cha.ch_som_vect().size_totale();
2021 }
2022
2023 return 0;
2024}
2025
2026int Champ_P1NC_implementation::nb_colonnes()
2027{
2028 if (sub_type(Champ_P1NC,le_champ()))
2029 {
2030 Champ_P1NC& cha =ref_cast(Champ_P1NC,le_champ());
2031 return cha.ch_som_vect().size_reelle();
2032 }
2033 else if (sub_type(Champ_Fonc_P1NC,le_champ()))
2034 {
2035 Champ_Fonc_P1NC& cha =ref_cast(Champ_Fonc_P1NC,le_champ());
2036 return cha.ch_som_vect().size_reelle();
2037 }
2038
2039 return 0;
2040}
2041
2042
2043void Champ_P1NC_implementation::dimensionner_Mat_Bloc_Morse_Sym(Matrice& matrice_tmp)
2044{
2045 // modif specific au Champ_P1NC_implementation
2046 const Matrice_Morse_Sym& la_matrice = get_MatP1NC2P1_L2();
2047 //
2048 int n1 = nb_colonnes_tot();
2049 int n2 = nb_colonnes();
2050 int iligne;
2051 const auto& tab1 = la_matrice.get_tab1();
2052 const auto& tab2 = la_matrice.get_tab2();
2053
2054 matrice_tmp.typer("Matrice_Bloc");
2055 Matrice_Bloc& matrice=ref_cast(Matrice_Bloc,matrice_tmp.valeur());
2056 matrice.dimensionner(1,2);
2057 matrice.get_bloc(0,0).typer("Matrice_Morse_Sym");
2058 matrice.get_bloc(0,1).typer("Matrice_Morse");
2059
2060 Matrice_Morse_Sym& MBrr = ref_cast(Matrice_Morse_Sym,matrice.get_bloc(0,0).valeur());
2061 Matrice_Morse& MBrv = ref_cast(Matrice_Morse,matrice.get_bloc(0,1).valeur());
2062 MBrr.dimensionner(n2,0);
2063 MBrv.dimensionner(n2,0);
2064
2065 auto& tab1RR = MBrr.get_set_tab1();
2066 auto& tab2RR = MBrr.get_set_tab2();
2067 auto& tab1RV = MBrv.get_set_tab1();
2068 auto& tab2RV = MBrv.get_set_tab2();
2069
2070 IntVect compteur_MBrr(n2);
2071 IntVect compteur_MBrv(n2);
2072 compteur_MBrr=0;
2073 compteur_MBrv=0;
2074
2075 // On parcours les lignes de la_matrice pour compter les elements
2076 // non nuls de chaque ligne
2077 int jcolonne;
2078 for (iligne=0; iligne<n2; iligne++)
2079 {
2080 for (auto k=tab1(iligne)-1; k<tab1(iligne+1)-1; k++)
2081 {
2082 jcolonne = tab2(k)-1;
2083 if (jcolonne < n2)
2084 {
2085 // l'element correspondant est dans la partie RR de la_matrice
2086 if ((jcolonne >= iligne) && (jcolonne < n2))
2087 {
2088 // l'element correspondant est situe au dessus de la diagonale de la_matrice
2089 compteur_MBrr(iligne)++;
2090 }
2091 }
2092 else
2093 {
2094 // l'element correspondant est dans la partie RV de la_matrice
2095 compteur_MBrv(iligne)++;
2096 }
2097 }
2098 }
2099
2100
2101 // On remplie tab1RR et tab1RV
2102 tab1RR(0)=1;
2103 tab1RV(0)=1;
2104 for(int i=0; i<n2; i++)
2105 {
2106 tab1RR(i+1)=compteur_MBrr(i)+tab1RR(i);
2107 tab1RV(i+1)=compteur_MBrv(i)+tab1RV(i);
2108 }
2109 // On dimensionne tab2RR et tab2RV
2110 MBrr.dimensionner(n2,tab1RR(n2)-1);
2111 MBrv.dimensionner(n2,n1-n2,tab1RV(n2)-1);
2112
2113 // On remplit tab2RR et tab2RV
2114 for (iligne=0; iligne<n2; iligne++)
2115 {
2116 auto compteurRR = tab1RR(iligne)-1;
2117 auto compteurRV = tab1RV(iligne)-1;
2118 for (auto k=tab1(iligne)-1; k<tab1(iligne+1)-1; k++)
2119 {
2120 jcolonne = tab2(k)-1;
2121 if (jcolonne < n2)
2122 {
2123 // l'element correspondant est dans la partie RR de la_matrice
2124 if ((jcolonne >= iligne) && (jcolonne < n2))
2125 {
2126 // l'element correspondant est situe au dessus de la diagonale de la_matrice
2127 tab2RR(compteurRR) = tab2(k);
2128 compteurRR++;
2129 }
2130 }
2131 else
2132 {
2133 // l'element correspondant est dans la partie RV de la_matrice
2134 tab2RV(compteurRV) = tab2(k)-n2;
2135 compteurRV++;
2136 }
2137 }
2138 }
2139}
2140
2141void Champ_P1NC_implementation::Mat_Morse_to_Mat_Bloc(Matrice& matrice_tmp)
2142{
2143 //
2144 const Matrice_Morse_Sym& la_matrice = get_MatP1NC2P1_L2();
2145 //
2146
2147 int n1 = nb_colonnes_tot();
2148 int n2 = nb_colonnes();
2149
2150 Matrice_Bloc& matrice= ref_cast(Matrice_Bloc,matrice_tmp.valeur());
2151 Matrice_Morse& MBrr = ref_cast(Matrice_Morse,matrice.get_bloc(0,0).valeur());
2152 Matrice_Morse& MBrv = ref_cast(Matrice_Morse,matrice.get_bloc(0,1).valeur());
2153
2154 auto& tab1RR = MBrr.get_set_tab1();
2155 auto& tab2RR = MBrr.get_set_tab2();
2156 auto& coeffRR = MBrr.get_set_coeff();
2157 auto& tab1RV = MBrv.get_set_tab1();
2158 auto& tab2RV = MBrv.get_set_tab2();
2159 auto& coeffRV = MBrv.get_set_coeff();
2160
2161 DoubleTab ligne_tmp(n1);
2162 for(int i=0; i<n2; i++)
2163 {
2164 // On recopie le premier bloc de la matrice dans un tableau :
2165 // ligne_tmp = 0;
2166 for (auto k=la_matrice.get_tab1()(i)-1; k<la_matrice.get_tab1()(i+1)-1; k++)
2167 ligne_tmp(la_matrice.get_tab2()(k) - 1) = la_matrice.get_coeff()(k);
2168
2169 // On complete la partie reelle de la matrice
2170 for (auto k=tab1RR(i)-1; k<tab1RR(i+1)-1; k++)
2171 coeffRR[k] = ligne_tmp(tab2RR[k] - 1);
2172
2173 // On complete la partie virtuelle
2174 for (auto k=tab1RV(i)-1; k<tab1RV(i+1)-1; k++)
2175 coeffRV[k] = ligne_tmp(n2 + tab2RV[k] - 1);
2176 }
2177
2178
2179 /* Cerr<<"Impression de la_matrice"<<finl;
2180 la_matrice.imprimer_formatte(Cerr);
2181 Cerr<<"Impression de MBrr"<<finl;
2182 MBrr.imprimer_formatte(Cerr);
2183 Cerr<<"Impression de MBrv"<<finl;
2184 MBrv.imprimer_formatte(Cerr);*/
2185
2186 /* Debog::verifier_Mat_faces("avant resolution systeme : la_matrice",la_matrice);
2187 Debog::verifier_Mat_faces("avant resolution systeme : MBrr",MBrr);*/
2188}
2189
Champ_Don_base::valeurs
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
Definition Champ_Don_base.h:49
Champ_Fonc_P1NC
classe Champ_Fonc_P1NC
Definition Champ_Fonc_P1NC.h:28
Champ_Fonc_P1NC::domaine_vef
const Domaine_VEF & domaine_vef() const override
Definition Champ_Fonc_P1NC.h:37
Champ_Inc_base::domaine
const Domaine & domaine() const
Definition Champ_Inc_base.cpp:648
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
Definition Champ_Inc_base.h:77
Champ_P1NC_implementation::fixer_nb_valeurs_nodales
int fixer_nb_valeurs_nodales(int)
Definition Champ_P1NC_implementation.cpp:1996
Champ_P1NC_implementation::valeur_a_sommet_compo
double valeur_a_sommet_compo(int num_som, int le_poly, int ncomp) const
Definition Champ_P1NC_implementation.cpp:1776
Champ_P1NC_implementation::valeur_aux_sommets
DoubleTab & valeur_aux_sommets(const Domaine &dom, DoubleTab &ch_som) const override
Definition Champ_P1NC_implementation.cpp:1659
Champ_P1NC_implementation::domaine_vef
virtual const Domaine_VEF & domaine_vef() const =0
Champ_P1NC_implementation::MatP1NC2P1_L2
Matrice_Morse_Sym MatP1NC2P1_L2
Definition Champ_P1NC_implementation.h:115
Champ_P1NC_implementation::nb_colonnes_tot
int nb_colonnes_tot()
Definition Champ_P1NC_implementation.cpp:2009
Champ_P1NC_implementation::ch_som_vect_
DoubleVect ch_som_vect_
Definition Champ_P1NC_implementation.h:127
Champ_P1NC_implementation::valeur_a_elem_compo
double valeur_a_elem_compo(const DoubleVect &position, int le_poly, int ncomp) const override
Definition Champ_P1NC_implementation.cpp:1346
Champ_P1NC_implementation::valeur_P1_H1
friend DoubleTab & valeur_P1_H1(const Champ_P1NC &, const Domaine &, DoubleTab &)
Definition Champ_P1NC_implementation.cpp:702
Champ_P1NC_implementation::imprime_P1NC
int imprime_P1NC(Sortie &, int) const
Definition Champ_P1NC_implementation.cpp:1974
Champ_P1NC_implementation::valeur_aux_elems_compo
DoubleVect & valeur_aux_elems_compo(const DoubleTab &positions, const IntVect &les_polys, DoubleVect &valeurs, int ncomp) const override
Definition Champ_P1NC_implementation.cpp:1492
Champ_P1NC_implementation::remplir_coord_noeuds
DoubleTab & remplir_coord_noeuds(DoubleTab &positions) const override
Definition Champ_P1NC_implementation.cpp:1941
Champ_P1NC_implementation::valeur_aux_elems_compo_smooth
DoubleVect & valeur_aux_elems_compo_smooth(const DoubleTab &positions, const IntVect &les_polys, DoubleVect &valeurs, int ncomp)
Definition Champ_P1NC_implementation.cpp:1604
Champ_P1NC_implementation::coord_barycentrique
double coord_barycentrique(const IntTab &sommet_poly, const DoubleTab &coord, double x, double y, int le_poly, int face) const
Definition Champ_P1NC_implementation.h:59
Champ_P1NC_implementation::fonction_forme_3D
double fonction_forme_3D(double x, double y, double z, int le_poly, int face, const IntTab &sommet_poly, const DoubleTab &coord) const
Definition Champ_P1NC_implementation.h:54
Champ_P1NC_implementation::filtrer_L2_deja_appele_
int filtrer_L2_deja_appele_
Definition Champ_P1NC_implementation.h:128
Champ_P1NC_implementation::Champ_P1NC_implementation
Champ_P1NC_implementation()
Definition Champ_P1NC_implementation.cpp:41
Champ_P1NC_implementation::dimensionner_Mat_Bloc_Morse_Sym
void dimensionner_Mat_Bloc_Morse_Sym(Matrice &matrice_tmp)
Definition Champ_P1NC_implementation.cpp:2043
Champ_P1NC_implementation::valeur_aux_sommets_compo
DoubleVect & valeur_aux_sommets_compo(const Domaine &dom, DoubleVect &ch_som, int ncomp) const override
Definition Champ_P1NC_implementation.cpp:1838
Champ_P1NC_implementation::ch_som_vect
DoubleVect & ch_som_vect()
Definition Champ_P1NC_implementation.h:100
Champ_P1NC_implementation::filtrer_H1
void filtrer_H1(DoubleTab &) const
Definition Champ_P1NC_implementation.cpp:881
Champ_P1NC_implementation::remplir_coord_noeuds_et_polys
int remplir_coord_noeuds_et_polys(DoubleTab &positions, IntVect &polys) const override
Definition Champ_P1NC_implementation.cpp:1958
Champ_P1NC_implementation::Mat_Morse_to_Mat_Bloc
void Mat_Morse_to_Mat_Bloc(Matrice &matrice_tmp)
Definition Champ_P1NC_implementation.cpp:2141
Champ_P1NC_implementation::filtrer_resu
void filtrer_resu(DoubleTab &) const
Definition Champ_P1NC_implementation.cpp:1110
Champ_P1NC_implementation::valeur_aux_elems_smooth
DoubleTab & valeur_aux_elems_smooth(const DoubleTab &positions, const IntVect &les_polys, DoubleTab &valeurs)
Definition Champ_P1NC_implementation.cpp:1537
Champ_P1NC_implementation::ch_som_
DoubleTab ch_som_
Definition Champ_P1NC_implementation.h:126
Champ_P1NC_implementation::fonction_forme_3D_v
KOKKOS_INLINE_FUNCTION double fonction_forme_3D_v(double x, double y, double z, int le_poly, int face, CIntTabView sommet_poly, CDoubleTabView coord) const
Definition Champ_P1NC_implementation.h:44
Champ_P1NC_implementation::valeur_aux_centres_de_gravite
DoubleTab & valeur_aux_centres_de_gravite(const Domaine &, DoubleTab &valeurs) const
Computes values at the centers of gravity for a P1NC field.
Definition Champ_P1NC_implementation.cpp:1393
Champ_P1NC_implementation::valeur_aux_elems
DoubleTab & valeur_aux_elems(const DoubleTab &positions, const IntVect &les_polys, DoubleTab &valeurs) const override
Definition Champ_P1NC_implementation.cpp:1430
Champ_P1NC_implementation::test
friend int test(Champ_P1NC &, const Domaine &)
Champ_P1NC_implementation::valeur_P1_L2
friend DoubleTab & valeur_P1_L2(Champ_P1NC &, const Domaine &)
Projection du champ P1NC "cha" sur l'espace des champs P1.
Definition Champ_P1NC_implementation.cpp:65
Champ_P1NC_implementation::MatP1NC2P1_H1
Matrice_Morse_Sym MatP1NC2P1_H1
Definition Champ_P1NC_implementation.h:120
Champ_P1NC_implementation::solveur_L2
SolveurSys solveur_L2
Definition Champ_P1NC_implementation.h:119
Champ_P1NC_implementation::ch_som
DoubleTab & ch_som()
Definition Champ_P1NC_implementation.h:99
Champ_P1NC_implementation::solveur_H1
SolveurSys solveur_H1
Definition Champ_P1NC_implementation.h:121
Champ_P1NC_implementation::get_MatP1NC2P1_L2
const Matrice_Morse_Sym & get_MatP1NC2P1_L2() const
Definition Champ_P1NC_implementation.h:102
Champ_P1NC_implementation::MatP1NC2P1_L2_Parallele
Matrice MatP1NC2P1_L2_Parallele
Definition Champ_P1NC_implementation.h:117
Champ_P1NC_implementation::nb_colonnes
int nb_colonnes()
Definition Champ_P1NC_implementation.cpp:2026
Champ_P1NC_implementation::fonction_forme_2D_v
KOKKOS_INLINE_FUNCTION double fonction_forme_2D_v(double x, double y, int le_poly, int face, CIntTabView sommet_poly, CDoubleTabView coord) const
Definition Champ_P1NC_implementation.h:39
Champ_P1NC_implementation::filtrer_L2
void filtrer_L2(DoubleTab &) const
Definition Champ_P1NC_implementation.cpp:931
Champ_P1NC_implementation::fonction_forme_2D
double fonction_forme_2D(double x, double y, int le_poly, int face, const IntTab &sommet_poly, const DoubleTab &coord) const
Definition Champ_P1NC_implementation.h:49
Champ_P1NC_implementation::valeur_a_elem
DoubleVect & valeur_a_elem(const DoubleVect &position, DoubleVect &val, int le_poly) const override
Definition Champ_P1NC_implementation.cpp:1317
Champ_P1NC
Definition Champ_P1NC.h:35
Champ_P1NC::valeur_a_sommet_compo
double valeur_a_sommet_compo(int num_som, int le_poly, int ncomp) const override
renvoi la compo eme corrdonne des valeurs a l'element le_poly au sommet sommet
Definition Champ_P1NC.h:77
Champ_P1NC::domaine_vef
const Domaine_VEF & domaine_vef() const override
Definition Champ_P1NC.h:66
Champ_Proto::valeurs
virtual DoubleTab & valeurs()=0
Champ_base
classe Champ_base Cette classe est la base de la hierarchie des champs.
Definition Champ_base.h:43
Champ_implementation::get_domaine_geom
const Domaine & get_domaine_geom() const
Definition Champ_implementation.cpp:80
Champ_implementation::le_champ
virtual Champ_base & le_champ()=0
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Conds_lim
classe Conds_lim Cette classe represente un vecteur de conditions aux limites.
Definition Conds_lim.h:32
Debog::verifier
static void verifier(const char *const msg, double)
Definition Debog.cpp:21
Dirichlet_homogene
Classe Dirichlet_homogene Cette classe est la classe de base de la hierarchie des conditions aux limi...
Definition Dirichlet_homogene.h:31
Dirichlet_paroi_fixe
classe Dirichlet_paroi_fixe Represente une paroi immobile dans une equation de type Navier_Stokes.
Definition Dirichlet_paroi_fixe.h:26
Domaine_32_64::nb_som_elem
int nb_som_elem() const
Renvoie le nombre de sommets des elements geometriques constituants le domaine.
Definition Domaine.h:474
Domaine_32_64::nb_elem_tot
int_t nb_elem_tot() const
Definition Domaine.h:132
Domaine_32_64::get_renum_som_perio
int_t get_renum_som_perio(int_t i) const
Definition Domaine.h:281
Domaine_32_64::les_elems
IntTab_t & les_elems()
Definition Domaine.h:129
Domaine_32_64::nb_faces_elem
int nb_faces_elem(int=0) const
Renvoie le nombre de face de type i des elements geometriques constituants le domaine.
Definition Domaine.h:484
Domaine_32_64::coord_sommets
const DoubleTab_t & coord_sommets() const
Definition Domaine.h:112
Domaine_32_64::creer_tableau_sommets
virtual void creer_tableau_sommets(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
Cree un tableau ayant une "ligne" par sommet du maillage.
Definition Domaine.cpp:1000
Domaine_32_64::nb_som_tot
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...
Definition Domaine.h:123
Domaine_32_64::nb_som
int_t nb_som() const
Renvoie le nombre de sommets du domaine.
Definition Domaine.h:121
Domaine_32_64::sommet_elem
int_t sommet_elem(int_t i, int j) const
Renvoie le numero (global) du j-ieme sommet du i-ieme element.
Definition Domaine.h:136
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:40
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Renvoie la i-ieme condition aux limites.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:54
Domaine_VF::nb_faces
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
Domaine_VF::volumes_entrelaces
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
Domaine_VF::nb_faces_tot
int nb_faces_tot() const
renvoie le nombre total de faces.
Definition Domaine_VF.h:481
Domaine_VF::face_normales
virtual double face_normales(int face, int comp) const
Definition Domaine_VF.h:47
Domaine_VF::xv
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
Domaine_VF::volumes
double volumes(int i) const
Definition Domaine_VF.h:113
Domaine_VF::face_sommets
int face_sommets(int i, int j) const
renvoie le numero du ieme sommet de la face num_face.
Definition Domaine_VF.h:583
Domaine_VF::nb_som_face
int nb_som_face() const
renvoie le nombre de sommets par face.
Definition Domaine_VF.h:494
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::oriente_normale
int oriente_normale(int f, int e) const
Definition Domaine_VF.h:194
Domaine_VF::premiere_face_int
int premiere_face_int() const
une face est interne ssi elle separe deux elements.
Definition Domaine_VF.h:463
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::nb_elem
int nb_elem() const
Definition Domaine_dis_base.h:49
Domaine_dis_base::nb_front_Cl
int nb_front_Cl() const
Definition Domaine_dis_base.h:53
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
EFichier
Fichier en lecture Cette classe est a la classe C++ ifstream ce que la classe Entree est a la.
Definition EFichier.h:29
Entree_Fichier_base::ouvrir
virtual int ouvrir(const char *name, IOS_OPEN_MODE mode=ios::in)
Definition Entree_Fichier_base.cpp:69
Equation_base::le_nom
const Nom & le_nom() const override
Renvoie le nom de l'equation.
Definition Equation_base.h:334
Equation_base::get_modele
virtual const RefObjU & get_modele(Type_modele type) const
Definition Equation_base.cpp:1894
Equation_base::domaine_Cl_dis
virtual Domaine_Cl_dis_base & domaine_Cl_dis()
Renvoie le domaine des conditions aux limite discretisee associee a l'equation.
Definition Equation_base.h:343
Field_base::le_nom
const Nom & le_nom() const override
Renvoie le nom du champ.
Definition Field_base.cpp:30
Field_base::nb_comp
virtual int nb_comp() const
Definition Field_base.h:56
Front_VF
class Front_VF
Definition Front_VF.h:36
Front_VF::nb_faces
int nb_faces() const
Definition Front_VF.h:53
Front_VF::num_premiere_face
int num_premiere_face() const
Definition Front_VF.h:63
Front_VF::nb_faces_tot
int nb_faces_tot() const
Definition Front_VF.h:58
Front_VF::num_face
int num_face(const int) const
Definition Front_VF.h:68
LecFicDiffuse
Cette classe implemente les operateurs et les methodes virtuelles de la classe EFichier de la facon s...
Definition LecFicDiffuse.h:32
LecFicDiffuse::ouvrir
int ouvrir(const char *name, IOS_OPEN_MODE mode=ios::in) override
Ouverture du fichier.
Definition LecFicDiffuse.cpp:61
Lec_Diffuse_base::good
int good() override
Definition Lec_Diffuse_base.cpp:100
Matrice_Bloc
Definition Matrice_Bloc.h:51
Matrice_Bloc::dimensionner
virtual void dimensionner(int N, int M)
Definition Matrice_Bloc.cpp:593
Matrice_Bloc::get_bloc
virtual const Matrice & get_bloc(int i, int j) const
Definition Matrice_Bloc.cpp:601
Matrice_Morse_Sym
Classe Matrice_Morse_Sym Represente une matrice M (creuse) symetrique stockee au format Morse.
Definition Matrice_Morse_Sym.h:34
Matrice_Morse_Sym::compacte
void compacte(int elim_coeff_nul=0)
Suppression des doublons on ordonne tab2;.
Definition Matrice_Morse_Sym.cpp:737
Matrice_Morse
Classe Matrice_Morse Represente une matrice M (creuse), non necessairement carree.
Definition Matrice_Morse.h:50
Matrice_Morse::get_set_tab2
auto & get_set_tab2()
Definition Matrice_Morse.h:103
Matrice_Morse::get_tab2
const auto & get_tab2() const
Definition Matrice_Morse.h:111
Matrice_Morse::dimensionner
void dimensionner(int n, _SIZE_ nnz)
Size the matrix with n lines and n columns and nnz zero-values coefficients.
Definition Matrice_Morse.cpp:305
Matrice_Morse::get_tab1
const auto & get_tab1() const
Definition Matrice_Morse.h:110
Matrice_Morse::get_set_coeff
auto & get_set_coeff()
Definition Matrice_Morse.h:108
Matrice_Morse::get_set_tab1
auto & get_set_tab1()
Definition Matrice_Morse.h:98
Matrice_Morse::get_coeff
const auto & get_coeff() const
Definition Matrice_Morse.h:112
Matrice_Morse::remplir
void remplir(const IntLists &, const DoubleLists &, const DoubleVect &)
Definition Matrice_Morse.cpp:1613
Matrice_Morse::nb_lignes
int nb_lignes() const override
Return local number of lines (=size on the current proc).
Definition Matrice_Morse.h:90
Matrice_Sym::set_est_definie
void set_est_definie(int)
Definition Matrice_Sym.cpp:25
Matrice
Classe Matrice Classe generique de la hierarchie des matrices.
Definition Matrice.h:34
Modele_turbulence_hyd_base
Classe Modele_turbulence_hyd_base Cette classe sert de base a la hierarchie des classes.
Definition Modele_turbulence_hyd_base.h:43
Modele_turbulence_hyd_base::utiliser_loi_paroi
bool utiliser_loi_paroi() const
Definition Modele_turbulence_hyd_base.h:52
Modele_turbulence_hyd_base::loi_paroi
const Turbulence_paroi_base & loi_paroi() const
Definition Modele_turbulence_hyd_base.h:50
MorEqn::equation
const Equation_base & equation() const
Renvoie la reference sur l'equation pointe par MorEqn::mon_equation.
Definition MorEqn.h:62
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
Periodique
classe Periodique Cette classe represente une condition aux limites periodique.
Definition Periodique.h:31
Periodique::face_associee
int face_associee(int i) const
Definition Periodique.h:35
Precond_base
Definition Precond_base.h:26
Process::mp_min
static double mp_min(double)
Definition Process.cpp:386
Process::me
static int me()
renvoie mon rang dans le groupe de communication courant.
Definition Process.cpp:125
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Solv_GCP
Definition Solv_GCP.h:26
SolveurSys
class SolveurSys Un SolveurSys represente n'importe qu'elle classe
Definition SolveurSys.h:32
SolveurSys::resoudre_systeme
int resoudre_systeme(const Matrice_Base &matrice, const DoubleVect &secmem, DoubleVect &solution)
Definition SolveurSys.cpp:48
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
TRUSTTab::ref
virtual void ref(const TRUSTTab &)
Definition TRUSTTab.tpp:308
TRUSTTab::view_wo
std::enable_if_t< is_default_exec_space< EXEC_SPACE >, View< _TYPE_, _SHAPE_ > > view_wo()
Definition TRUSTTab.h:276
TRUSTTab::nb_dim
int nb_dim() const
Definition TRUSTTab.h:199
TRUSTTab::resize
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469
TRUSTTab::dimension_tot
_SIZE_ dimension_tot(int) const override
Definition TRUSTTab.tpp:160
TRUSTTab::view_ro
std::enable_if_t< is_default_exec_space< EXEC_SPACE >, ConstView< _TYPE_, _SHAPE_ > > view_ro() const
Definition TRUSTTab.h:261
TRUSTTab::view_rw
std::enable_if_t< is_default_exec_space< EXEC_SPACE >, View< _TYPE_, _SHAPE_ > > view_rw()
Definition TRUSTTab.h:291
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133
TRUSTVect::size
_SIZE_ size() const
Definition TRUSTVect.tpp:45
TRUSTVect::size_totale
_SIZE_ size_totale() const
Definition TRUSTVect.tpp:61
TRUSTVect::line_size
int line_size() const
Definition TRUSTVect.tpp:67
TRUSTVect::size_reelle
_SIZE_ size_reelle() const
Definition TRUSTVect.tpp:27
TRUSTVect::resize
void resize(_SIZE_, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTVect.tpp:91
TRUSTVect::echange_espace_virtuel
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")
Definition TRUSTVect.tpp:282
TRUST_Ref_Objet_U::valeur
const Objet_U & valeur() const
Definition TRUST_Ref.h:134
Turbulence_paroi_base
Classe Turbulence_paroi_base Classe de base pour la hierarchie des classes representant les modeles.
Definition Turbulence_paroi_base.h:38
Turbulence_paroi_base::Cisaillement_paroi
const DoubleTab & Cisaillement_paroi() const
Definition Turbulence_paroi_base.h:126