TrioCFD 1.9.9_beta
TrioCFD documentation
Loading...
Searching...
No Matches
Quadri_VEF.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 <Quadri_VEF.h>
17#include <Domaine.h>
18#include <Domaine_VEF.h>
19
20Implemente_instanciable_sans_constructeur(Quadri_VEF,"Quadri_VEF",Elem_VEF_base);
21
22// printOn and readOn
23
24Sortie& Quadri_VEF::printOn(Sortie& s ) const
25{
26 return s << que_suis_je() << finl;
27}
28
29Entree& Quadri_VEF::readOn(Entree& s )
30{
31 return s ;
32}
33
34/*! @brief KEL_(0,fa7),KEL_(1,fa7) are the local indices of the 2 faces surrounding facet with local index fa7.
35 *
36 * The local index of fa7 is that of the node that carries it.
37 *
38 */
39Quadri_VEF::Quadri_VEF()
40{
41 int tmp[3][4]=
42 {
43 {0, 1, 0, 2},
44 {1, 2, 3, 3},
45 {0, 1, 2, 3}
46 };
47 KEL_.resize(3,4);
48 for (int i=0; i<3; i++)
49 for (int j=0; j<4; j++)
50 KEL_(i,j)=tmp[i][j];
51}
52
53/*! @brief Fills the face_normales array in Domaine_VEF.
54 *
55 */
56void Quadri_VEF::creer_face_normales(DoubleTab& Face_normales,
57 const IntTab& Face_sommets,
58 const IntTab& Face_voisins,
59 const IntTab& elem_faces,
60 const Domaine& domaine_geom) const
61{
62 const DoubleTab& les_coords = domaine_geom.coord_sommets();
63 int nb_face_tot = Face_normales.dimension_tot(0);
64 for (int num_Face=0; num_Face<nb_face_tot; num_Face++)
65 {
66 double x1, y1;
67 double nx, ny;
68 double x1g = 0, y1g = 0;
69 double x2g = 0, y2g = 0;
70 double grx, gry, psc;
71 int sign = 1, i;
72 int n0 = Face_sommets(num_Face, 0);
73 int n1 = Face_sommets(num_Face, 1);
74 x1 = les_coords(n0, 0) - les_coords(n1, 0);
75 y1 = les_coords(n0, 1) - les_coords(n1, 1);
76 nx = -y1;
77 ny = x1;
78 int elem1 = Face_voisins(num_Face, 0);
79 int elem2 = Face_voisins(num_Face, 1);
80
81 //Orient the normal towards the element with the higher index
82 //for that, first check whether we are on the boundary
83 if (elem2 != -1)
84 {
85 //orient from the center of gravity
86 //compute the center of gravity of each element
87 for (i = 0; i < 4; i++)
88 {
89 x1g += les_coords(Face_sommets(elem_faces(elem1, i), 0), 0);
90 x1g += les_coords(Face_sommets(elem_faces(elem1, i), 1), 0);
91 y1g += les_coords(Face_sommets(elem_faces(elem1, i), 0), 1);
92 y1g += les_coords(Face_sommets(elem_faces(elem1, i), 1), 1);
93 x2g += les_coords(Face_sommets(elem_faces(elem2, i), 0), 0);
94 x2g += les_coords(Face_sommets(elem_faces(elem2, i), 1), 0);
95 y2g += les_coords(Face_sommets(elem_faces(elem2, i), 0), 1);
96 y2g += les_coords(Face_sommets(elem_faces(elem2, i), 1), 1);
97 }
98
99 grx = (x2g - x1g) * 0.125;
100 gry = (y2g - y1g) * 0.125;
101
102 //check the sign of the dot product
103 psc = grx * nx + gry * ny;
104 if (psc < 0)
105 {
106 if (elem1 < elem2)
107 sign = -1;
108 else if (elem2 < elem1)
109 sign = -1;
110 }
111 }
112 else
113 {
114 //orient from the center of gravity and the midpoint of the
115 //current face
116
117 for (i = 0; i < 4; i++)
118 {
119 x1g += les_coords(Face_sommets(elem_faces(elem1, i), 0), 0);
120 x1g += les_coords(Face_sommets(elem_faces(elem1, i), 1), 0);
121 y1g += les_coords(Face_sommets(elem_faces(elem1, i), 0), 1);
122 y1g += les_coords(Face_sommets(elem_faces(elem1, i), 1), 1);
123 }
124 // Cerr << "xg et yg de Face_normales: " << x1g << " " << y1g << finl;
125
126 x2g = les_coords(n0, 0) + les_coords(n1, 0);
127 y2g = les_coords(n0, 1) + les_coords(n1, 1);
128 grx = x2g * 0.5 - x1g * 0.125;
129 gry = y2g * 0.5 - y1g * 0.125;
130
131 // Cerr << "grx et gry : " << grx << " " << gry << finl;
132 //check the sign of the dot product
133 psc = grx * nx + gry * ny;
134 if (psc < 0)
135 sign = -1;
136 }
137 Face_normales(num_Face, 0) = sign * nx;
138 Face_normales(num_Face, 1) = sign * ny;
139 }
140}
141
142/*! @brief Computes the facet normals for standard elements.
143 *
144 */
145void Quadri_VEF::creer_facette_normales(const Domaine_VEF& dom_VEF,
146 const IntVect& rang_elem_non_std) const
147{
148 const Domaine& domaine_geom = dom_VEF.domaine();
149 auto& facette_normales = const_cast<Domaine_VEF&>(dom_VEF).facette_normales();
150 const DoubleTab& les_coords = domaine_geom.coord_sommets();
151 const IntTab& les_Polys = domaine_geom.les_elems();
152 int nb_elem_tot = domaine_geom.nb_elem_tot();
153
154 int i, fa7;
155 int i0=-1,i1=-1;
156 int num_som[4];
157 double x[4][2];
158 double xg[2];
159 double xj0[2];
160 double u[2];
161 double v[2];
162 double psc;
163
164 if (facette_normales.dimension(0) != nb_elem_tot)
165 facette_normales.resize(nb_elem_tot,nb_facette(),2);
166
167 for(i=0; i<nb_elem_tot; i++)
168 {
169 if (rang_elem_non_std(i)==-1)
170 {
171 num_som[0]=les_Polys(i,0);
172 num_som[1]=les_Polys(i,1);
173 num_som[2]=les_Polys(i,2);
174 num_som[3]=les_Polys(i,3);
175 // Cerr << "numero de sommet dans Quadri_VEF : "
176 // << num_som[0] << " " << num_som[1] << " "
177 // << num_som[2] << " " << num_som[3] << finl;
178 x[0][0]=les_coords(num_som[0],0);
179 x[0][1]=les_coords(num_som[0],1);
180 x[1][0]=les_coords(num_som[1],0);
181 x[1][1]=les_coords(num_som[1],1);
182 x[2][0]=les_coords(num_som[2],0);
183 x[2][1]=les_coords(num_som[2],1);
184 x[3][0]=les_coords(num_som[3],0);
185 x[3][1]=les_coords(num_som[3],1);
186 xg[0]=(x[0][0]+x[1][0]+x[2][0]+x[3][0])*0.25;
187 xg[1]=(x[0][1]+x[1][1]+x[2][1]+x[3][1])*0.25;
188
189 for (fa7=0; fa7<4; fa7++)
190 {
191 // facet fa7 of a standard element has vertices
192 // fa7 and G with coordinates xg:
193 u[0]= x[fa7][0]-xg[0];
194 u[1]= x[fa7][1]-xg[1];
195 v[0]= -u[1];
196 v[1]= u[0];
197
198 // Orient normals:
199 switch (fa7)
200 {
201 case 0:
202 i0=1;
203 i1=2;
204 break;
205 case 1:
206 i0=3;
207 i1=0;
208 break;
209 case 2:
210 i0=3;
211 i1=0;
212 break;
213 case 3:
214 i0=2;
215 i1=1;
216 break;
217 default:
218 i0=-1;
219 i1=-1;
220 assert(0);
221 exit();
222 break;
223 }
224
225 // i0 = KEL_(0,fa7);
226 // i1 = KEL_(1,fa7);
227
228 xj0[0]= x[i0][0]-x[i1][0];
229 xj0[1]= x[i0][1]-x[i1][1];
230
231 psc=xj0[0]*v[0] + xj0[1]*v[1];
232 // Cerr << "facette_normales : " << facette_normales << finl;
233 //Cerr << "numero d'element : " << i << finl;
234 //Cerr << "numero de fa7 : " << fa7 << finl;
235 if (psc < 0)
236 {
237 facette_normales(i,fa7,0) = -v[0];
238 facette_normales(i,fa7,1) = -v[1];
239 }
240 else
241 {
242 facette_normales(i,fa7,0) = v[0];
243 facette_normales(i,fa7,1) = v[1];
244 }
245 //Cerr << "facette_normales : " << facette_normales << finl;
246 }
247 }
248 }
249
250 // Cerr << "facette_normales internes (Quadri_VEF::creer_normales_facettes) : " << facette_normales << finl;
251}
252
253/*! @brief Fills the normales_facettes_Cl array in Domaine_Cl_VEF for facet fa7 of element num_elem.
254 *
255 */
256void Quadri_VEF::creer_normales_facettes_Cl(DoubleTab& normales_facettes_Cl,
257 int fa7,
258 int num_elem,const DoubleTab& x,
259 const DoubleVect& xg, const Domaine& domaine_geom) const
260{
261 double u[2];
262 double v[2];
263 double xj0[2];
264 double psc;
265
266 u[0]= x(fa7,0)-xg[0];
267 u[1]= x(fa7,1)-xg[1];
268 v[0]= -u[1];
269 v[1]= u[0];
270
271 int i0=-1;
272 int i1=-1;
273
274 switch (fa7)
275 {
276 case 0:
277 i0=1;
278 i1=2;
279 break;
280 case 1:
281 i0=3;
282 i1=0;
283 break;
284 case 2:
285 i0=3;
286 i1=0;
287 break;
288 case 3:
289 i0=2;
290 i1=1;
291 break;
292 default:
293 i0=-1;
294 i1=-1;
295 assert(0);
296 exit();
297 break;
298 }
299
300 // Orient normals:
301 xj0[0]= x(i0,0)-x(i1,0);
302 xj0[1]= x(i0,1)-x(i1,1);
303
304 psc=xj0[0]*v[0] + xj0[1]*v[1];
305 if (psc < 0)
306 {
307 normales_facettes_Cl(num_elem,fa7,0) = -v[0];
308 normales_facettes_Cl(num_elem,fa7,1) = -v[1];
309 }
310 else
311 {
312 normales_facettes_Cl(num_elem,fa7,0) = v[0];
313 normales_facettes_Cl(num_elem,fa7,1) = v[1];
314 }
315 // Cerr << "normales_facettes_Cl " << normales_facettes_Cl << finl;
316}
317
318/*! @brief Modifies the interlaced volumes for face j of a non-standard element.
319 *
320 */
321void Quadri_VEF::modif_volumes_entrelaces(int j,int elem,
322 const Domaine_VEF& le_dom_VEF,
323 DoubleVect& volumes_entrelaces_Cl,
324 int type_cl) const
325{
326
327 Cerr << "Quadri_VEF::modif_volumes_entrelaces() does nothing for now" << finl;
328 // double surf_mod;
329 // const DoubleVect& volumes_entrelaces = le_dom_VEF.volumes_entrelaces();
330 // const IntTab& elem_faces = le_dom_VEF.elem_faces();
331 // Cerr << "elem " << elem << finl;
332 // Cerr << "type_cl " << type_cl << finl;
333
334 // switch(type_cl) {
335
336 // // no Dirichlet face: impossible
337 // case 0:
338 // {
339 // Cerr << "Quadri_VEF::modif_volumes_entrelaces() type 0 impossible!\n";
340 // break;
341 // }
342
343 // case 1: // one Dirichlet face: Face 3
344 // {
345 // DoubleTab coord(3,2);
346 // for(int i=0;i<3;i++)
347 // for(int k=0;k<2;k++)
348 // coord(i,k)=le_dom_VEF.domaine().coord(le_dom_VEF.domaine().sommet_elem(elem,i),k);
349 // volumes_entrelaces_Cl[elem_faces(elem,0)] = std::fabs((coord(0,0)-coord(2,0))*((coord(2,1)+coord(3,1))/2-coord(2,2)) - ((coord(2,0)+coord(3,0))/2-coord(2,0))*(coord(0,1)-coord(2,2)))/2;
350 // volumes_entrelaces_Cl[elem_faces(elem,2)] = std::fabs((coord(1,0)-coord(3,0))*((coord(2,1)+coord(3,1))/2-coord(3,2)) - ((coord(2,0)+coord(3,0))/2-coord(3,0))*(coord(1,1)-coord(3,2)))/2;
351 // volumes_entrelaces_Cl[elem_faces(elem,1)] = std::fabs((coord(1,0)-coord(0,0))*((coord(2,1)+coord(3,1))/2-coord(0,2)) - ((coord(2,0)+coord(3,0))/2-coord(0,0))*(coord(1,1)-coord(0,2)))/2;
352 // break;
353 // }
354
355 // case 3: // one Dirichlet face: Face 2
356 // {
357 // DoubleTab coord(3,2);
358 // for(int i=0;i<3;i++)
359 // for(int k=0;k<2;k++)
360 // coord(i,k)=le_dom_VEF.domaine().coord(le_dom_VEF.domaine().sommet_elem(elem,i),k);
361 // volumes_entrelaces_Cl[elem_faces(elem,1)] = std::fabs((coord(2,0)-coord(3,0))*((coord(1,1)+coord(3,1))/2-coord(3,2)) - ((coord(1,0)+coord(3,0))/2-coord(3,0))*(coord(2,1)-coord(3,2)))/2;
362 // volumes_entrelaces_Cl[elem_faces(elem,3)] = std::fabs((coord(0,0)-coord(1,0))*((coord(1,1)+coord(3,1))/2-coord(1,2)) - ((coord(1,0)+coord(3,0))/2-coord(1,0))*(coord(0,1)-coord(1,2)))/2;
363 // volumes_entrelaces_Cl[elem_faces(elem,0)] = std::fabs((coord(0,0)-coord(2,0))*((coord(1,1)+coord(3,1))/2-coord(2,2)) - ((coord(1,0)+coord(3,0))/2-coord(2,0))*(coord(0,1)-coord(2,2)))/2;
364 // break;
365 // }
366
367 // case 9: // one Dirichlet face: Face 1
368 // {
369 // DoubleTab coord(3,2);
370 // for(int i=0;i<3;i++)
371 // for(int k=0;k<2;k++)
372 // coord(i,k)=le_dom_VEF.domaine().coord(le_dom_VEF.domaine().sommet_elem(elem,i),k);
373 // volumes_entrelaces_Cl[elem_faces(elem,0)] = std::fabs((coord(2,0)-coord(0,0))*((coord(0,1)+coord(1,1))/2-coord(0,2)) - ((coord(0,0)+coord(1,0))/2-coord(0,0))*(coord(2,1)-coord(0,2)))/2;
374 // volumes_entrelaces_Cl[elem_faces(elem,2)] = std::fabs((coord(3,0)-coord(1,0))*((coord(0,1)+coord(1,1))/2-coord(1,2)) - ((coord(0,0)+coord(1,0))/2-coord(1,0))*(coord(3,1)-coord(1,2)))/2;
375 // volumes_entrelaces_Cl[elem_faces(elem,3)] = std::fabs((coord(3,0)-coord(2,0))*((coord(0,1)+coord(1,1))/2-coord(2,2)) - ((coord(0,0)+coord(1,0))/2-coord(2,0))*(coord(3,1)-coord(2,2)))/2;
376 // break;
377 // }
378
379 // case 27: // one Dirichlet face: Face 0
380 // {
381
382 // DoubleTab coord(3,2);
383 // for(int i=0;i<3;i++) {
384 // for(int k=0;k<2;k++) {
385 // Cerr << "le_dom_VEF.domaine().sommet_elem(elem,i) " << le_dom_VEF.domaine().sommet_elem(elem,i) << finl;
386 // coord(i,k)=le_dom_VEF.domaine().coord(le_dom_VEF.domaine().sommet_elem(elem,i),k);
387 // Cerr << "coord(i,k) " << coord(i,k) << finl;
388 // }
389 // }
390 // Cerr << " elem_faces(elem,1) " << elem_faces(elem,1) << finl;
391 // Cerr << " elem_faces(elem,3) " << elem_faces(elem,3) << finl;
392 // Cerr << " elem_faces(elem,2) " << elem_faces(elem,2) << finl;
393
394 // Cerr << " coord(0,0) " << coord(0,0) << " coord(0,1) " << coord(0,1) << finl;
395 // Cerr << " coord(1,0) " << coord(1,0) << " coord(1,1) " << coord(1,1) << finl;
396 // Cerr << " coord(2,0) " << coord(2,0) << " coord(2,1) " << coord(2,1) << finl;
397
398 // volumes_entrelaces_Cl[elem_faces(elem,1)] = std::fabs((coord(3,0)-coord(2,0))*((coord(0,1)+coord(2,1))/2-coord(2,2)) - ((coord(0,0)+coord(2,0))/2-coord(2,0))*(coord(3,1)-coord(2,2)))/2;
399 // volumes_entrelaces_Cl[elem_faces(elem,3)] = std::fabs((coord(1,0)-coord(0,0))*((coord(0,1)+coord(2,1))/2-coord(0,2)) - ((coord(0,0)+coord(2,0))/2-coord(0,0))*(coord(1,1)-coord(0,2)))/2;
400 // volumes_entrelaces_Cl[elem_faces(elem,2)] = std::fabs((coord(3,0)-coord(1,0))*((coord(0,1)+coord(2,1))/2-coord(1,2)) - ((coord(0,0)+coord(2,0))/2-coord(1,0))*(coord(3,1)-coord(1,2)))/2;
401
402 // break;
403 // }
404
405 // case 4: // two Dirichlet faces: Faces 2,3
406 // {
407 // Cerr << " elem_faces(elem,0) " << elem_faces(elem,0) << finl;
408 // Cerr << " elem_faces(elem,1) " << elem_faces(elem,1) << finl;
409 // volumes_entrelaces_Cl[elem_faces(elem,0)] +=volumes_entrelaces[elem_faces(elem,3)];
410 // volumes_entrelaces_Cl[elem_faces(elem,1)] +=volumes_entrelaces[elem_faces(elem,2)];
411 // break;
412 // }
413
414 // case 28: // two Dirichlet faces: Faces 0,3
415 // {
416 // volumes_entrelaces_Cl[elem_faces(elem,1)] +=volumes_entrelaces[elem_faces(elem,0)];
417 // volumes_entrelaces_Cl[elem_faces(elem,2)] +=volumes_entrelaces[elem_faces(elem,3)];
418 // break;
419 // }
420
421 // case 12: // two Dirichlet faces: Faces 1,2
422 // {
423 // volumes_entrelaces_Cl[elem_faces(elem,0)] +=volumes_entrelaces[elem_faces(elem,1)];
424 // volumes_entrelaces_Cl[elem_faces(elem,3)] +=volumes_entrelaces[elem_faces(elem,2)];
425 // break;
426 // }
427
428 // case 36: // two Dirichlet faces: Faces 0,1
429 // {
430 // volumes_entrelaces_Cl[elem_faces(elem,3)] +=volumes_entrelaces[elem_faces(elem,0)];
431 // volumes_entrelaces_Cl[elem_faces(elem,2)] +=volumes_entrelaces[elem_faces(elem,1)];
432 // break;
433 // }
434
435 // case 10: // two Dirichlet faces: Faces 1,3
436 // {
437 // volumes_entrelaces_Cl[elem_faces(elem,0)] +=volumes_entrelaces[elem_faces(elem,0)];
438 // volumes_entrelaces_Cl[elem_faces(elem,2)] +=volumes_entrelaces[elem_faces(elem,2)];
439 // break;
440 // }
441
442 // case 30: // two Dirichlet faces: Faces 0,2
443 // {
444 // volumes_entrelaces_Cl[elem_faces(elem,1)] +=volumes_entrelaces[elem_faces(elem,1)];
445 // volumes_entrelaces_Cl[elem_faces(elem,3)] +=volumes_entrelaces[elem_faces(elem,3)];
446 // break;
447 // }
448
449 // case 13: //three Dirichlet faces: Faces 3,2,1
450 // {
451 // volumes_entrelaces_Cl[elem_faces(elem,0)] +=volumes_entrelaces[elem_faces(elem,3)]+volumes_entrelaces[elem_faces(elem,2)]
452 // +volumes_entrelaces[elem_faces(elem,1)];
453 // break;
454 // }
455
456 // case 31: //three Dirichlet faces: Faces 0,3,2
457 // {
458 // volumes_entrelaces_Cl[elem_faces(elem,1)] +=volumes_entrelaces[elem_faces(elem,3)]+volumes_entrelaces[elem_faces(elem,2)]
459 // +volumes_entrelaces[elem_faces(elem,0)];
460 // break;
461 // }
462
463 // case 37: //three Dirichlet faces: Faces 1,0,3
464 // {
465 // volumes_entrelaces_Cl[elem_faces(elem,2)] +=volumes_entrelaces[elem_faces(elem,3)]+volumes_entrelaces[elem_faces(elem,1)]
466 // +volumes_entrelaces[elem_faces(elem,0)];
467 // break;
468 // }
469
470 // case 39: //three Dirichlet faces: Faces 2,1,0
471 // {
472 // volumes_entrelaces_Cl[elem_faces(elem,3)] +=volumes_entrelaces[elem_faces(elem,2)]+volumes_entrelaces[elem_faces(elem,1)]
473 // +volumes_entrelaces[elem_faces(elem,0)];
474 // break;
475 // }
476
477 // default :
478 // {
479 // Cerr << "\n type inconnu : " << type_cl ;
480 // exit();
481 // }
482
483 // } // end switch
484
485}
486
487/*! @brief Modifies the interlaced volumes for joint face j of a non-standard element.
488 *
489 */
490void Quadri_VEF::modif_volumes_entrelaces_faces_joints(int j,int elem,
491 const Domaine_VEF& le_dom_VEF,
492 DoubleVect& volumes_entrelaces_Cl,
493 int type_cl) const
494{
495 // const DoubleVect& volumes_entrelaces = le_dom_VEF.volumes_entrelaces();
496 // const IntTab& elem_faces = le_dom_VEF.elem_faces();
497
498 // switch(type_cl) {
499
500 // // no Dirichlet face: impossible
501 // case 0:
502 // {
503 // Cerr << "Quadri_VEF::modif_volumes_entrelaces() type 0 impossible!\n";
504 // break;
505 // }
506
507
508 // } // end switch
509
510}
511
512/*! @brief
513 *
514 */
515void Quadri_VEF::calcul_vc(const ArrOfInt& Face,ArrOfDouble& vc,
516 const ArrOfDouble& vs,const DoubleTab& vsom,
517 const Champ_Inc_base& vitesse,int type_cl, const DoubleVect& porosite_face) const
518{
519 //Cerr << " DANS Quadri_VEF::calcul_vc , type_cl = " << type_cl << finl;
520 //Cerr << "vs " << vs << " et vsom " << vsom << " et vitesse " << vitesse << finl;
521
522 // switch(type_cl) {
523 // case 0: // no Dirichlet face
524 // {
525 vc[0] = vs[0]*0.25;
526 vc[1] = vs[1]*0.25;
527 // break;
528 // }
529
530 // case 1: // one Dirichlet face: Face 3
531 // {
532 // vc[0] = vitesse.valeurs()(Face[3],0);
533 // vc[1] = vitesse.valeurs()(Face[3],1);
534 // break;
535 // }
536
537 // case 3: // one Dirichlet face: Face 2
538 // {
539 // vc[0] = vitesse.valeurs()(Face[2],0);
540 // vc[1] = vitesse.valeurs()(Face[2],1);
541 // break;
542 // }
543
544 // case 9: // one Dirichlet face: Face 1
545 // {
546 // vc[0] = vitesse.valeurs()(Face[1],0);
547 // vc[1] = vitesse.valeurs()(Face[1],1);
548 // break;
549 // }
550
551 // case 27: // one Dirichlet face: Face 0
552 // {
553 // vc[0] = vitesse.valeurs()(Face[0],0);
554 // vc[1] = vitesse.valeurs()(Face[0],1);
555 // break;
556 // }
557
558 // case 4: // two Dirichlet faces: Faces 2,3
559 // {
560 // vc[0]= vsom(3,0);
561 // vc[1]= vsom(3,1);
562 // break;
563 // }
564
565 // case 28: // two Dirichlet faces: Faces 0,3
566 // {
567 // vc[0]= vsom(2,0);
568 // vc[1]= vsom(2,1);
569 // break;
570 // }
571
572 // case 12: // two Dirichlet faces: Faces 1,2
573 // {
574 // vc[0]= vsom(1,0);
575 // vc[1]= vsom(1,1);
576 // break;
577 // }
578
579 // case 36: // two Dirichlet faces: Faces 0,1
580 // {
581 // vc[0]= vsom(0,0);
582 // vc[1]= vsom(0,1);
583 // break;
584 // }
585
586
587 // case 10: // two Dirichlet faces: Faces 1,3
588 // {
589 // vc[0] = vs[0]*0.25;
590 // vc[1] = vs[1]*0.25;
591 // break;
592 // }
593
594 // case 30: // two Dirichlet faces: Faces 0,2
595 // {
596 // vc[0] = vs[0]*0.25;
597 // vc[1] = vs[1]*0.25;
598 // break;
599 // }
600
601 // case 13: //three Dirichlet faces: Faces 3,2,1
602 // {
603 // vc[0]= vitesse.valeurs()(Face[2],0);
604 // vc[1]= vitesse.valeurs()(Face[2],1);
605 // break;
606 // }
607
608 // case 31: //three Dirichlet faces: Faces 0,3,2
609 // {
610 // vc[0]= vitesse.valeurs()(Face[3],0);
611 // vc[1]= vitesse.valeurs()(Face[3],1);
612 // break;
613 // }
614
615 // case 37: //three Dirichlet faces: Faces 1,0,3
616 // {
617 // vc[0]= vitesse.valeurs()(Face[0],0);
618 // vc[1]= vitesse.valeurs()(Face[0],1);
619 // break;
620 // }
621
622 // case 39: //three Dirichlet faces: Faces 2,1,0
623 // {
624 // vc[0]= vitesse.valeurs()(Face[1],0);
625 // vc[1]= vitesse.valeurs()(Face[1],1);
626 // break;
627 // }
628
629 // default :
630 // {
631 // Cerr << "\n type inconnu : " << type_cl ;
632 // exit();
633 // }
634
635 // } // end switch
636
637}
638
639/*! @brief Computes the coordinates xg of the centre of a non-standard element.
640 *
641 * Also computes idirichlet = number of Dirichlet faces of the element.
642 * If idirichlet=2, n1 is the index of the node coinciding with G.
643 *
644 */
645void Quadri_VEF::calcul_xg(DoubleVect& xg, const DoubleTab& x,
646 const int type_elem_Cl,int& idirichlet,int& n1,int& ,int& ) const
647{
648 int dim=xg.size();
649 // switch(type_elem_Cl) {
650
651 // case 0: // no Dirichlet face: it has 4 facets
652 // // point G is the barycentre of the triangle vertices
653 // {
654 for (int j=0; j<dim; j++)
655 xg[j]=(x(0,j)+x(1,j)+x(2,j)+x(3,j))*0.25;
656 idirichlet=0;
657 // break;
658 // }
659
660 // case 1: // one Dirichlet face: Face 3
661 // // point G is the center of face 3
662 // {
663 // for (int j=0; j<dim; j++)
664 // xg[j]=(x(2,j)+x(3,j))*0.5;
665 // idirichlet=1;
666 // break;
667 // }
668
669 // case 3: // one Dirichlet face: Face 2
670 // // point G is the center of face 2
671 // {
672 // for (int j=0; j<dim; j++)
673 // xg[j]=(x(1,j)+x(3,j))*0.5;
674 // idirichlet=1;
675 // break;
676 // }
677
678 // case 9: // one Dirichlet face: Face 1
679 // // point G is the center of face 1
680 // {
681 // for (int j=0; j<dim; j++)
682 // xg[j]=(x(0,j)+x(1,j))*0.5;
683 // idirichlet=1;
684 // break;
685 // }
686
687 // case 27: // one Dirichlet face: Face 0
688 // // point G is the center of face 0
689 // {
690 // for (int j=0; j<dim; j++)
691 // xg[j]=(x(0,j)+x(2,j))*0.5;
692 // idirichlet=1;
693 // break;
694 // }
695
696 // case 4: // two Dirichlet faces: Faces 2,3
697 // //point G is the vertex common to the two Dirichlet faces
698 // {
699 // for (int j=0; j<dim; j++)
700 // xg[j]=x(3,j);
701 // idirichlet=2;
702 // break;
703 // }
704
705 // case 28: // two Dirichlet faces: Faces 0,3
706 // //point G is the vertex common to the two Dirichlet faces
707 // {
708 // for (int j=0; j<dim; j++)
709 // xg[j]=x(2,j);
710 // idirichlet=2;
711 // break;
712 // }
713
714 // case 12: // two Dirichlet faces: Faces 1,2
715 // //point G is the vertex common to the two Dirichlet faces
716 // {
717 // for (int j=0; j<dim; j++)
718 // xg[j]=x(1,j);
719 // idirichlet=2;
720 // break;
721 // }
722
723 // case 36: // two Dirichlet faces: Faces 0,1
724 // //point G is the vertex common to the two Dirichlet faces
725 // {
726 // for (int j=0; j<dim; j++)
727 // xg[j]=x(0,j);
728 // idirichlet=2;
729 // break;
730 // }
731
732 // case 10: // two Dirichlet faces: Faces 1,3
733 // // keep the same control volumes as for internal faces
734 // {
735 // for (int j=0; j<dim; j++)
736 // xg[j]=(x(0,j)+x(1,j)+x(2,j)+x(3,j))*0.25;
737 // idirichlet=0;
738 // break;
739 // }
740
741 // case 30: // two Dirichlet faces: Faces 0,2
742 // // keep the same control volumes as for internal faces
743 // {
744 // for (int j=0; j<dim; j++)
745 // xg[j]=(x(0,j)+x(1,j)+x(2,j)+x(3,j))*0.25;
746 // idirichlet=0;
747 // break;
748 // }
749
750 // case 13: //three Dirichlet faces: Faces 3,2,1
751 // // point G is the center of the Dirichlet face opposite to the non-Dirichlet face
752 // {
753 // for (int j=0; j<dim; j++)
754 // xg[j]=(x(1,j)+x(3,j))*0.5;
755 // idirichlet=3;
756 // break;
757 // }
758
759 // case 31: //three Dirichlet faces: Faces 0,3,2
760 // // point G is the center of the Dirichlet face opposite to the non-Dirichlet face
761 // {
762 // for (int j=0; j<dim; j++)
763 // xg[j]=(x(2,j)+x(3,j))*0.5;
764 // idirichlet=3;
765 // break;
766 // }
767
768 // case 37: //three Dirichlet faces: Faces 1,0,3
769 // // point G is the center of the Dirichlet face opposite to the non-Dirichlet face
770 // {
771 // for (int j=0; j<dim; j++)
772 // xg[j]=(x(0,j)+x(2,j))*0.5;
773 // idirichlet=3;
774 // break;
775 // }
776
777 // case 39: //three Dirichlet faces: Faces 2,1,0
778 // // point G is the center of the Dirichlet face opposite to the non-Dirichlet face
779 // {
780 // for (int j=0; j<dim; j++)
781 // xg[j]=(x(0,j)+x(1,j))*0.5;
782 // idirichlet=3;
783 // break;
784 // }
785
786 // } // end switch
787
788}
789
790
791/*! @brief Modifies normales_facettes_Cl when idirichlet=3.
792 *
793 * idirichlet = number of Dirichlet faces of the element.
794 * If idirichlet=3, n1 is the index of the node coinciding with G.
795 *
796 */
797void Quadri_VEF::modif_normales_facettes_Cl(DoubleTab& normales_facettes_Cl,
798 int fa7,int num_elem,
799 int idirichlet,int n1,int ,int ) const
800{
801 // Cerr << "Quadri_VEF::modif_normales_facettes_Cl, fa7 " << fa7 << "num_elem " << num_elem << "idirichlet " << idirichlet << "n1 " << n1 << "normales_facettes_Cl " << normales_facettes_Cl << finl;
802 switch (idirichlet)
803 {
804 case 0:
805 break;
806 case 1:
807 break;
808 case 2:
809 break;
810 case 3:
811 {
812 // normales_facettes_Cl(num_elem,n1,0) = 0;
813 // normales_facettes_Cl(num_elem,n1,1) = 0;
814 break;
815 }
816 }
817}
Champ_Inc_base
Class Champ_Inc_base.
Definition Champ_Inc_base.h:53
Domaine_32_64::nb_elem_tot
int_t nb_elem_tot() const
Definition Domaine.h:132
Domaine_32_64::les_elems
IntTab_t & les_elems()
Definition Domaine.h:129
Domaine_32_64::coord_sommets
const DoubleTab_t & coord_sommets() const
Definition Domaine.h:112
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:53
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Elem_VEF_base
Definition Elem_VEF_base.h:27
Elem_VEF_base::KEL_
IntTab KEL_
Definition Elem_VEF_base.h:41
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Objet_U::que_suis_je
const Nom & que_suis_je() const
Returns the string identifying the class.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Reads an Objet_U from an input stream. Virtual method to override.
Definition Objet_U.cpp:289
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Writes the object to an output stream. Virtual method to override.
Definition Objet_U.cpp:278
Process::exit
static void exit(int exit_code=-1)
Exit routine for TRUST within a Kokkos region.
Definition Process.cpp:466
Quadri_VEF
Definition Quadri_VEF.h:24
Quadri_VEF::creer_normales_facettes_Cl
void creer_normales_facettes_Cl(DoubleTab &, int, int, const DoubleTab &, const DoubleVect &, const Domaine &) const override
Fills the normales_facettes_Cl array in Domaine_Cl_VEF for facet fa7 of element num_elem.
Definition Quadri_VEF.cpp:256
Quadri_VEF::nb_facette
int nb_facette() const override
Definition Quadri_VEF.h:30
Quadri_VEF::calcul_vc
void calcul_vc(const ArrOfInt &, ArrOfDouble &, const ArrOfDouble &, const DoubleTab &, const Champ_Inc_base &, int, const DoubleVect &) const override
Definition Quadri_VEF.cpp:515
Quadri_VEF::modif_volumes_entrelaces_faces_joints
void modif_volumes_entrelaces_faces_joints(int, int, const Domaine_VEF &, DoubleVect &, int) const override
Modifies the interlaced volumes for joint face j of a non-standard element.
Definition Quadri_VEF.cpp:490
Quadri_VEF::creer_face_normales
void creer_face_normales(DoubleTab &, const IntTab &, const IntTab &, const IntTab &, const Domaine &) const override
Fills the face_normales array in Domaine_VEF.
Definition Quadri_VEF.cpp:56
Quadri_VEF::creer_facette_normales
void creer_facette_normales(const Domaine_VEF &, const IntVect &) const override
Computes the facet normals for standard elements.
Definition Quadri_VEF.cpp:145
Quadri_VEF::modif_volumes_entrelaces
void modif_volumes_entrelaces(int, int, const Domaine_VEF &, DoubleVect &, int) const override
Modifies the interlaced volumes for face j of a non-standard element.
Definition Quadri_VEF.cpp:321
Quadri_VEF::calcul_xg
void calcul_xg(DoubleVect &, const DoubleTab &, const int, int &, int &, int &, int &) const override
Computes the coordinates xg of the centre of a non-standard element.
Definition Quadri_VEF.cpp:645
Quadri_VEF::modif_normales_facettes_Cl
void modif_normales_facettes_Cl(DoubleTab &, int, int, int, int, int, int) const override
Modifies normales_facettes_Cl when idirichlet=3.
Definition Quadri_VEF.cpp:797
Quadri_VEF::Quadri_VEF
Quadri_VEF()
KEL_(0,fa7),KEL_(1,fa7) are the local indices of the 2 faces surrounding facet with local index fa7.
Definition Quadri_VEF.cpp:39
Sortie
Base class for output streams.
Definition Sortie.h:52
TRUSTTab::dimension_tot
_SIZE_ dimension_tot(int) const override
Definition TRUSTTab.tpp:160
TRUSTVect::size
_SIZE_ size() const
Definition TRUSTVect.tpp:45