TrioCFD 1.9.9_beta
TrioCFD documentation
Loading...
Searching...
No Matches
Assembleur_P_VEFPre1B_tools.cpp
1/****************************************************************************
2* Copyright (c) 2025, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
12* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
13*
14*****************************************************************************/
15
16#include <Assembleur_P_VEFPreP1B.h>
17#include <Matrice_Bloc_Sym.h>
18#include <Domaine.h>
19#include <Dirichlet.h>
20#include <Dirichlet_homogene.h>
21#include <Periodique.h>
22#include <Neumann_sortie_libre.h>
23#include <Symetrie.h>
24#include <EcrFicCollecteBin.h>
25#include <TRUSTLists.h>
26#include <Navier_Stokes_std.h>
27#include <Op_Div_VEFP1B_Elem.h>
28#include <Op_Grad_VEF_P1B_Face.h>
29#include <Milieu_base.h>
30#include <Scatter.h>
31#include <communications.h>
32
33static int face_associee=-1;
34
35static ArrOfDouble gradi(3);
36static ArrOfDouble gradj(3);
37static inline void
38projette(ArrOfDouble& grad, int face, const DoubleTab& normales)
39{
40 double psc=0, norm=0;
41 int dimension=Objet_U::dimension, comp;
42 for(comp=0; comp<dimension; comp++)
43 {
44 psc+=grad[comp]*normales(face,comp);
45 norm+=normales(face,comp)*normales(face,comp);
46 }
47 // psc/=norm; // Fixed bug: Arithmetic exception
48 if (std::fabs(norm)>=DMINFLOAT) psc/=norm;
49 for(comp=0; comp<dimension; comp++)
50 {
51 grad[comp]-=psc*normales(face,comp);
52 }
53 psc=0;
54 // for(comp=0; comp<dimension; comp++)
55 // {
56 // psc+=gradi(comp)*normales(face,comp);
57 // }
58 // assert(psc < 1.e-10);
59}
60//
61
62// returns the first non-Dirichlet face
63// 2 if Perio
64// 3 if Neumann
65// 4 if Symmetry
66// 1 otherwise
67// and face_associee=-1 except if periodic (face_associee=associated face)
68static inline int okface(int& ind_face, int& face, const Cond_lim& la_cl)
69{
70 face_associee=-1;
71 int ok=1;
72 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
73 int nb_faces_bord_tot = le_bord.nb_faces_tot();
74 do
75 {
76 face=le_bord.num_face(ind_face);
77 if ((sub_type(Dirichlet, la_cl.valeur()))
78 || (sub_type(Dirichlet_homogene, la_cl.valeur())))
79 {
80 ok=0;
81 }
82 else if (sub_type(Periodique,la_cl.valeur()))
83 {
84 //periodicity
85 const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
86 face_associee=le_bord.num_face(la_cl_perio.face_associee(ind_face));
87 ok=2;
88 }
89 else if (sub_type(Neumann_sortie_libre, la_cl.valeur()))
90 {
91 //sortie_libre
92 ok=3;
93 }
94 else if (sub_type(Symetrie, la_cl.valeur()))
95 {
96 //symmetry
97 ok=4;
98 }
99 }
100 while ( ( (ok==0) || ((ok==2)&&(face_associee<face)) ) && (++ind_face<nb_faces_bord_tot) );
101 if (ind_face==nb_faces_bord_tot) ok=-1;
102 return ok;
103}
104
105inline int verifier_complet(const Assembleur_P_VEFPreP1B& ass,
106 const Matrice_Bloc_Sym& matrice,
107 const Domaine_VEF& domaine_VEF)
108{
109 const Navier_Stokes_std& eqn=ref_cast(Navier_Stokes_std, ass.equation());
110 const Operateur_Div& opdiv=eqn.operateur_divergence();
111 const Op_Div_VEFP1B_Elem& div=ref_cast(Op_Div_VEFP1B_Elem,
112 opdiv.valeur());
113 //div.verifier();
114 const Operateur_Grad& opgrad=eqn.operateur_gradient();
115 const Op_Grad_VEF_P1B_Face& grad=ref_cast(Op_Grad_VEF_P1B_Face,
116 opgrad.valeur());
117 //grad.verifier();
118 const Solveur_Masse_base& solvm=eqn.solv_masse();
119 const DoubleTab& pression=eqn.pression().valeurs();
120 DoubleTab tab(pression);
121 DoubleTab resu(tab), resu2(tab);
122
123 // Compute an arbitrary pressure field
124 exemple_champ_non_homogene(domaine_VEF, tab);
125
126 DoubleTab gradP(eqn.inconnue().valeurs());
127 grad.calculer(tab, gradP);
128 solvm.appliquer(gradP);
129 div.calculer(gradP, resu);
130 matrice.multvect(tab, resu2);
131 return 0;
132}
133int verifier( const Assembleur_P_VEFPreP1B& ass,
134 const Matrice_Bloc_Sym& matrice,
135 const Domaine_VEF& domaine_VEF,
136 const DoubleTab& inverse_quantitee_entrelacee)
137{
138 if (domaine_VEF.get_alphaE()+domaine_VEF.get_alphaS()+domaine_VEF.get_alphaA()!=Objet_U::dimension)
139 {
140 Cerr << "Assembleur_P_VEFPreP1B::verifier is not designed to verify your discretization." << finl;
141 Process::exit();
142 }
143 const Navier_Stokes_std& eqn=ref_cast(Navier_Stokes_std, ass.equation());
144 const Operateur_Div& opdiv=eqn.operateur_divergence();
145 const Op_Div_VEFP1B_Elem& div=ref_cast(Op_Div_VEFP1B_Elem,
146 opdiv.valeur());
147 //div.verifier();
148 const Operateur_Grad& opgrad=eqn.operateur_gradient();
149 const Op_Grad_VEF_P1B_Face& grad=ref_cast(Op_Grad_VEF_P1B_Face,
150 opgrad.valeur());
151 //grad.verifier();
152 // const Solveur_Masse_base& solvm=eqn.solv_masse();
153 const DoubleTab& pression=eqn.pression().valeurs();
154 int ko=0;
155 ko=verifier_complet(ass, matrice, domaine_VEF);
156 DoubleTab pre(pression);
157 DoubleTab resu(pre), resu2(pre), erreur(pre);
158 DoubleTab gradP(eqn.inconnue().valeurs());
159 const Domaine& domaine=domaine_VEF.domaine();
160 int nb_elem=domaine.nb_elem();
161 int nb_elem_tot=domaine.nb_elem_tot();
162 int nb_som=domaine.nb_som();
163 int nb_som_tot=domaine.nb_som_tot();
164 int nb_aretes=domaine.nb_aretes();
165 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
166 for (int proc=0; proc<Process::nproc(); proc++)
167 {
168 int n = pre.dimension_tot(0);
169 // Processor proc broadcasts its value of n to all processes
170 envoyer_broadcast(n, proc);
171
172 for(int i=0; i<n; i++)
173 {
174 //Cerr << "[" << Process::me() << "] Checking row " << i << " of the matrix." << finl;
175 pre=0;
176 if (Process::me()==proc)
177 {
178 if (0<=i && i<nb_elem)
179 {
180 Cerr << "Checking real element " << i << " of processor " << proc;
181 pre(i)=1;
182 }
183 else if (nb_elem_tot<=i && i<nb_elem_tot+nb_som)
184 {
185 int sommet=i-nb_elem_tot;
186 Cerr << "Checking real vertex ";
187 if (domaine.get_renum_som_perio(sommet)!=sommet) Cerr << "periodic ";
188 Cerr << i-nb_elem_tot << " of processor " << proc;
189 pre(i)=1;
190 }
191 else if (nb_elem_tot+nb_som_tot<=i && i<nb_elem_tot+nb_som_tot+nb_aretes && ok_arete[i-nb_elem_tot-nb_som_tot])
192 {
193 int arete=i-nb_elem_tot-nb_som_tot;
194 const ArrOfInt& renum_arete_perio=domaine_VEF.get_renum_arete_perio();
195 if (renum_arete_perio[arete]==arete)
196 {
197 Cerr << "Checking real, non-redundant, non-periodic edge " << arete << " of processor " << proc;
198 pre(i)=1;
199 }
200 }
201 Cerr << " ";
202 //Cerr << finl;
203 }
204 double verifie=mp_max_vect(pre);
205 //verifie=1;pre(i)=1; // Test everything
206 if (verifie)
207 {
208 pre.echange_espace_virtuel();
209 /* The debug check is useless here: sequential and parallel results are not comparable.
210 Nom ch;
211 ch="pre for row ";
212 ch+=(Nom)i+" of processor ";
213 ch+=(Nom)proc+" =";
214 Debog::verifier(ch,pre);
215 */
216 // Compute via Div(Grad(P))
217 grad.calculer(pre, gradP);
218 {
219 int nbf=inverse_quantitee_entrelacee.size();
220 int d = Objet_U::dimension;
221 for (int face=0; face<nbf; face++)
222 for(int k=0; k<d; k++)
223 gradP(face,k)*=inverse_quantitee_entrelacee(face,k);
224 }
225 //solvm.appliquer(gradP);
226 div.calculer(gradP, resu);
227 // Compute via -Lap(P)
228 matrice.multvect(pre, resu2);
229 // The error should be zero
230 erreur=resu2;
231 erreur+=resu;
232 resu*=-1;
233 // Case where the diagonal is *2, correct it:
234 if (est_egal(2*resu(i),resu2(i))) erreur(i)=0;
235 // Optimization: combine 3 mp_norme_vect into 1 collective call
236 double erreur_carre = local_carre_norme_vect(erreur);
237 double resu2_carre = local_carre_norme_vect(resu2);
238 double resu_carre = local_carre_norme_vect(resu);
239 Process::mp_sum_for_each(erreur_carre, resu2_carre, resu_carre);
240 double erreur_absolue = sqrt(erreur_carre);
241 double erreur_relative = erreur_absolue / (sqrt(resu2_carre) + sqrt(resu_carre) + DMINFLOAT);
242 double app=mp_prodscal(resu,pre);
243 if(erreur_absolue>1.e-12 && erreur_relative>1.e-6)
244 {
245 Cerr << "[" << Process::me() << "] FAIL at row " << i << " for proc " << proc << " (AP,P)= " << app << " error= " << erreur_absolue << finl;
246 ko=1;
247 Cerr << "[" << Process::me() << "] pre= ";
248 pre.ecrit(Cerr);
249 Cerr << "[" << Process::me() << "] Div(gradP) = ";
250 resu.ecrit(Cerr);
251 Cerr << "[" << Process::me() << "] Lap(P) = ";
252 resu2.ecrit(Cerr);
253 Cerr << "[" << Process::me() << "] error = ";
254 erreur.ecrit(Cerr);
255 Cerr << "[" << Process::me() << "] invqtentrelacee = ";
256 inverse_quantitee_entrelacee.ecrit(Cerr);
257 }
258 else
259 {
260 Cerr << "[" << Process::me() << "] OK at row " << i << " for proc " << proc << " (AP,P)= " << app << " error= " << erreur_absolue << finl;
261 }
262 }
263 }
264 }
265 if (ko)
266 {
267 Cerr << "[" << Process::me() << "] Pressure matrix:" << finl;
268 matrice.imprimer_formatte(Cerr);
269 Cerr << "Failure in the assembler's verifier method. Check the FAIL entries above." << finl;
270 Process::exit();
271 }
272 return 1;
273}
274//
275// sort the sommets array in ascending order and
276// faces_op1 and faces_op2 accordingly.
277//
278static inline void sort( ArrOfInt& sommets, ArrOfInt& faces_op1, ArrOfInt& faces_op2)
279{
280 int sz=sommets.size_array();
281 if(sommets[sz-1]==-1) sz--;
282 int i,j;
283 for(i=0; i<sz; i++)
284 for(j=i; j<sz; j++)
285 if(sommets[i]>sommets[j])
286 {
287 int tmp=sommets[i];
288 sommets[i]=sommets[j];
289 sommets[j]=tmp;
290 tmp=faces_op1[i];
291 faces_op1[i]=faces_op1[j];
292 faces_op1[j]=tmp;
293 tmp=faces_op2[i];
294 faces_op2[i]=faces_op2[j];
295 faces_op2[j]=tmp;
296 }
297}
298static inline int chercher_arete(const Domaine_VEF& domaine_VEF,
299 int elem, int somi, int somj,
300 const IntTab& elem_aretes,
301 const IntTab& aretes_som)
302{
303 const ArrOfInt& renum_arete_perio=domaine_VEF.get_renum_arete_perio();
304 const Domaine& dom=domaine_VEF.domaine();
305 if(somi>somj)
306 {
307 int k=somi;
308 somi=somj;
309 somj=k;
310 }
311 for(int i_arete=0; i_arete<6; i_arete++)
312 {
313 int arete=elem_aretes(elem, i_arete);
314 int som1=dom.get_renum_som_perio(aretes_som(arete,0));
315 int som2=dom.get_renum_som_perio(aretes_som(arete,1));
316 somi=dom.get_renum_som_perio(somi);
317 somj=dom.get_renum_som_perio(somj);
318 if( (somi==som1)
319 && (somj==som2) )
320 return renum_arete_perio[arete];
321 if( (somi==som2)
322 && (somj==som1) )
323 return renum_arete_perio[arete];
324 }
325 return -1;
326}
327static inline void swap (int& i, int& j)
328{
329 int k=i;
330 i=j;
331 j=k;
332}
333
334//
335// Fill sommets, faces_op1 and faces_op2,
336// where sommets contains the vertices of face and the vertices of elem1 and elem2
337// that are not in face. face_op1(i) is the index of the face opposite to sommets(i)
338// in elem1. face_op2(i) is ... in elem2. (if elem2=-1, then face_op2=-1)
339// The first 'dimension' entries are the vertices of face
340// the last one is in elem2
341static inline void remplir_sommets(const Domaine_VEF& domaine_VEF,
342 int face, int elem1, int elem2,
343 ArrOfInt& sommets,
344 ArrOfInt& faces_op1,
345 ArrOfInt& faces_op2)
346{
347 int dplusun=Objet_U::dimension+1;
348 const IntTab& elem_som = domaine_VEF.domaine().les_elems();
349 const IntTab& face_som = domaine_VEF.face_sommets();
350 const IntTab& elem_faces = domaine_VEF.elem_faces();
351 const Domaine& dom=domaine_VEF.domaine();
352 for(int i=0; i<Objet_U::dimension; i++)
353 sommets[i]=dom.get_renum_som_perio(face_som(face,i));
354 if(elem1!=-1)
355 {
356 int ok=0;
357 for(int i=0; i<dplusun; i++)
358 if( (elem_faces(elem1,i)==face) ||
359 (elem_faces(elem1,i)==face_associee) )
360 {
361 sommets[Objet_U::dimension]=
362 dom.get_renum_som_perio(elem_som(elem1, i));
363 faces_op1[Objet_U::dimension]=face;
364 faces_op2[Objet_U::dimension]=-1;
365 ok=1;
366 }
367 else
368 {
369 int j=dom.get_renum_som_perio(elem_som(elem1, i));
370 for(int k=0; k<Objet_U::dimension; k++)
371 if(j==sommets[k])
372 faces_op1[k]=elem_faces(elem1, i);
373 }
374 if (ok!=1)
375 {
376 Cerr << "The discretization used has a P1 component" << finl;
377 Cerr << "which is not available to deal your mesh." << finl;
378 Cerr << "The mesh with this discretization must contain only ";
379 Cerr << (Objet_U::dimension==2?"triangles":"tetraedras") << "." << finl;
380 Process::exit();
381 }
382 }
383 else
384 {
385 Cerr << "Not handled: elem1==-1 case" << finl;
386 Process::exit();
387 }
388 if(elem2!=-1)
389 {
390 //int ok=0;
391 for(int i=0; i<dplusun; i++)
392 if( (elem_faces(elem2,i)==face)||
393 (elem_faces(elem2,i)==face_associee) )
394 {
395 sommets[dplusun]=dom.get_renum_som_perio(elem_som(elem2, i));
396 faces_op2[dplusun]=face;
397 faces_op1[dplusun]=-1;
398 //ok=1;
399 }
400 else
401 {
402 int j=dom.get_renum_som_perio(elem_som(elem2, i));
403 for(int k=0; k<Objet_U::dimension; k++)
404 if(j==sommets[k])
405 faces_op2[k]=elem_faces(elem2, i);
406 }
407 // assert(ok==1) is commented out because ok=1 is commented out
408 }
409 else
410 {
411 sommets[dplusun]=-1;
412 faces_op2[dplusun]=-1;
413 faces_op1[dplusun]=-1;
414 }
415}
416
417// compute the gradient at the face separating elem1 and elem2
418// of the shape function associated with vertex s
419//
420static void calculer_grad(const IntTab& face_voisins,
421 int elem1, int elem2,
422 const ArrOfDouble& coef_som,
423 int s, int fop1, int fop2,
424 const DoubleTab& normales,
425 ArrOfDouble& grad)
426{
427 int dimension=Objet_U::dimension;
428 double signe=1;
429 if(fop1!=-1)
430 {
431 if(elem1!=face_voisins(fop1,0))
432 signe=-1;
433 signe*=coef_som[elem1];
434 for(int comp=0; comp<dimension; comp++)
435 grad[comp]=signe*normales(fop1,comp);
436 }
437 else
438 grad=0;
439 if((elem2!=-1)&&(fop2!=-1))
440 {
441 signe=1;
442 if(elem2!=face_voisins(fop2,0))
443 signe=-1;
444 signe*=coef_som[elem2];
445 for(int comp=0; comp<dimension; comp++)
446 grad[comp]+=signe*normales(fop2,comp);
447 }
448}
449
450// compute the gradient at the face separating elem1 and elem2
451// of the shape function associated with vertex s
452//
453static void calculer_grad_arete(int face,
454 const IntTab& face_voisins,
455 int i, int j,
456 int elem1, int elem2,
457 int fop1, int fop2,
458 int fop3, int fop4,
459 const DoubleTab& normales,
460 ArrOfDouble& grad)
461{
462 assert(face_voisins(face,0)==elem1);
463 int signe1=1,signe2=1,signe3=1,signe4=1;
464 if((!(fop1==-1) && !(face_voisins(fop1,0)==elem1)))
465 signe1=-1;
466 if(!(fop3==-1) && !(face_voisins(fop3,0)==elem1))
467 signe3=-1;
468 if(elem2!=-1)
469 {
470 if((!(fop2==-1) && !(face_voisins(fop2,0)==elem2)))
471 signe2=-1;
472 if(!(fop4==-1) && !( face_voisins(fop4,0)==elem2))
473 signe4=-1;
474 }
475 if(j<3) // an edge of the face
476 {
477 if(elem2!=-1)
478 {
479 for(int comp=0; comp<3; comp++)
480 grad[comp]=-2./15.*(signe1*normales(fop1,comp)
481 +signe2*normales(fop2,comp)
482 +signe3*normales(fop3,comp)
483 +signe4*normales(fop4,comp));
484 }
485 else
486 {
487 for(int comp=0; comp<3; comp++)
488 grad[comp]=-2./15.*(signe1*normales(fop1,comp)
489 +signe3*normales(fop3,comp)
490 +normales(face,comp));
491 }
492
493 }
494 else if (j==3) // an edge of elem1 but not of a face
495 {
496 for(int comp=0; comp<3; comp++)
497 grad[comp]=1./15.*(signe1*normales(fop1,comp)
498 +signe3*normales(fop3,comp));
499 }
500 else // an edge of elem2 but not of a face
501 {
502 assert(j==4);
503 for(int comp=0; comp<3; comp++)
504 grad[comp]=1./15.*(signe2*normales(fop2,comp)
505 +signe4*normales(fop4,comp));
506 }
507}
508
509static double dotproduct_array_fois_inverse_quantitee_entrelacee(const ArrOfDouble& grad1,const ArrOfDouble& grad2,const DoubleTab& inverse_quantitee_entrelacee, int face )
510{
511 double dot=0;
512 int size=inverse_quantitee_entrelacee.dimension(1);
513 for (int i=0; i<size; i++) dot+=grad1[i]*grad2[i]*inverse_quantitee_entrelacee(face,i);
514 return dot;
515 //return dotproduct_array(gradi,gradj)*inverse_quantitee_entrelacee(face,0);
516}
517
518static void contribuer_matriceP0P1(int elem1, int elem2, const ArrOfInt& sommets,
519 IntLists& voisins, DoubleLists& coeffs)
520{
521 int dimension=Objet_U::dimension,
522 dplusdeux=dimension+2;
523
524 for(int i=0; i<dplusdeux; i++)
525 {
526 int si=sommets[i];
527 if (si<0) break;
528 int rang1=voisins[elem1].rang(si);
529 if(rang1==-1)
530 {
531 voisins[elem1].add(si);
532 coeffs[elem1].add(0);
533 }
534 if (elem2!=-1)
535 {
536 int rang2=voisins[elem2].rang(si);
537 if(rang2==-1)
538 {
539 voisins[elem2].add(si);
540 coeffs[elem2].add(0);
541 }
542 }
543 }
544}
545inline void range(int& i, int& n, int& j, int& m, Matrice_Morse& ARR, Matrice_Morse& ARV, Matrice_Morse& AVR, Matrice_Morse& AVV, double coeff)
546{
547 if(i<n)
548 if(j<m)
549 ARR(i,j)+=coeff;
550 else
551 ARV(i,j-m)+=coeff;
552 else if(j<m)
553 AVR(i-n,j)+=coeff;
554 else
555 AVV(i-n,j-m)+=coeff;
556}
557static void update_matriceP0P1(const Domaine_VEF& domaine_VEF,
558 const DoubleTab& inverse_quantitee_entrelacee,
559 int face, int elem1, int elem2,
560 ArrOfInt& sommets, ArrOfInt& faces_op1,
561 ArrOfInt& faces_op2, const ArrOfDouble& coef_som,
562 Matrice_Morse& ARR, Matrice_Morse& ARV,
563 Matrice_Morse& AVR, Matrice_Morse& AVV)
564{
565 const DoubleTab& normales = domaine_VEF.face_normales();
566 const IntTab& face_voisins = domaine_VEF.face_voisins();
567 assert(elem1==face_voisins(face, 0));
568 assert(elem2==face_voisins(face, 1));
569
570 int dimension=Objet_U::dimension,
571 dplusdeux=dimension+2;
572 double psc;
573 //double coeff_som=1./(dimension)/(dimension+1);
574
575 int nb_elem=domaine_VEF.nb_elem();
576 int nb_som=domaine_VEF.nb_som();
577 for(int i=0; i<dplusdeux; i++)
578 {
579 int si=sommets[i];
580 if (si<0) break;
581 calculer_grad(face_voisins, elem1, elem2, coef_som,si, faces_op1[i],
582 faces_op2[i], normales, gradi);
583 for(int k=0; k<dimension; k++)
584 gradj[k]=normales(face,k);
585 psc=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradj,
586 inverse_quantitee_entrelacee,face);
587 range(elem1,nb_elem,si,nb_som,ARR,ARV,AVR,AVV,psc);
588 if (elem2!=-1)
589 range(elem2,nb_elem,si,nb_som,ARR,ARV,AVR,AVV,-psc);
590 }
591}
592
593static void contribuer_matriceP1P1(int elem1, int elem2, const ArrOfInt& sommets, IntLists& voisins, DoubleLists& coeffs)
594
595{
596 int dimension=Objet_U::dimension,
597 dplusdeux=dimension+2;
598
599 // Skip vertices -1 which
600 // are at the end of the sommets array:
601 while (sommets[dplusdeux-1]==-1)
602 dplusdeux--;
603
604 for(int i=0; i<dplusdeux; i++)
605 {
606 int si=sommets[i];
607 for(int j=i+1; j<dplusdeux; j++)
608 {
609 int sj=sommets[j];
610 int rang=voisins[si].rang(sj);
611 if (sj>si)
612 {
613 if(rang==-1)
614 {
615 voisins[si].add(sj);
616 coeffs[si].add(0);
617 }
618 }
619 }
620 }
621}
622
623static void update_matriceP1P1(const Domaine_VEF& domaine_VEF,
624 const DoubleTab& inverse_quantitee_entrelacee,
625 int face, int elem1, int elem2,
626 ArrOfInt& sommets, ArrOfInt& faces_op1,
627 ArrOfInt& faces_op2, const ArrOfDouble& coef_som,
628 Matrice_Morse& ARR, Matrice_Morse& ARV,
629 Matrice_Morse& AVR, Matrice_Morse& AVV)
630{
631 const DoubleTab& normales = domaine_VEF.face_normales();
632 const IntTab& face_voisins = domaine_VEF.face_voisins();
633
634 int dimension=Objet_U::dimension,
635 dplusdeux=dimension+2;
636 double psc;
637 //double coeff_som=1./(dimension)/(dimension+1);
638 //coeff_som*=coeff_som;
639 int nb_som_tot=domaine_VEF.nb_som();
640 int i,j;
641
642 // Skip vertices -1 which
643 // are at the end of the sommets array:
644 while (sommets[dplusdeux-1]==-1)
645 dplusdeux--;
646
647 for(i=0; i<dplusdeux; i++)
648 {
649 int si=sommets[i];
650 calculer_grad(face_voisins, elem1, elem2, coef_som,si, faces_op1[i],
651 faces_op2[i], normales, gradi);
652 if(si<nb_som_tot)
653 ARR(si,si)+=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradi,
654 inverse_quantitee_entrelacee,face);
655 for(j=i+1; j<dplusdeux; j++)
656 {
657 int sj=sommets[j];
658 calculer_grad(face_voisins, elem1, elem2,coef_som, sj, faces_op1[j],
659 faces_op2[j], normales, gradj);
660 psc=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradj,inverse_quantitee_entrelacee,face);
661 //assert(sj>si);
662 if(si<nb_som_tot)
663 if(sj<nb_som_tot)
664 ARR(si,sj)+=psc;
665 else
666 ARV(si,sj-nb_som_tot)+=psc;
667 else if(sj<nb_som_tot)
668 AVR(si-nb_som_tot,sj)+=psc;
669 else
670 AVV(si-nb_som_tot,sj-nb_som_tot)+=psc;
671 }
672 }
673}
674
675static void contribuer_matricePaPa(const Domaine_VEF& domaine_VEF,
676 int elem1, int elem2,
677 const ArrOfInt& sommets, IntLists& voisins, DoubleLists& coeffs)
678
679{
680 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
681 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
682 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
683 int jmax=5;
684 if(elem2==-1) jmax=4;
685 for(int i=0; i<3; i++)
686 {
687 int si=sommets[i];
688 for(int j=i+1; j<jmax; j++)
689 {
690 int sj=sommets[j];
691 int arete1;
692 if(j<4)
693 arete1 = chercher_arete(domaine_VEF,elem1, si, sj,
694 elem_aretes, aretes_som);
695 else
696 arete1 = chercher_arete(domaine_VEF,elem2, si, sj,
697 elem_aretes, aretes_som);
698 assert(arete1!=-1);
699 if(ok_arete[arete1])
700 {
701 int jj=j;
702 for(int k=i; k<3; k++)
703 {
704 int sk=sommets[k];
705 for(int l=jj+1; l<jmax; l++)
706 {
707 int sl=sommets[l];
708 int arete2;
709 if(l<4)
710 arete2 = chercher_arete(domaine_VEF,elem1, sl, sk,
711 elem_aretes,
712 aretes_som);
713 else
714 arete2 = chercher_arete(domaine_VEF,elem2, sl, sk,
715 elem_aretes,
716 aretes_som);
717 assert(arete2!=-1);
718 if(ok_arete[arete2])
719 {
720 int tmp=arete1;
721 if(arete1>arete2) swap(arete1, arete2);
722 int rang=voisins[arete1].rang(arete2);
723 if(rang==-1)
724 {
725 voisins[arete1].add(arete2);
726 coeffs[arete1].add(0);
727
728 }
729 arete1=tmp;
730 }
731 }
732 jj=k+1;
733 }
734 }
735 }
736 }
737}
738static void update_matricePaPa(const Domaine_VEF& domaine_VEF,
739 const DoubleTab& inverse_quantitee_entrelacee,
740 int face, int elem1, int elem2,
741 ArrOfInt& sommets, ArrOfInt& faces_op1,
742 ArrOfInt& faces_op2,
743 Matrice_Morse& ARR, Matrice_Morse& ARV,
744 Matrice_Morse& AVR, Matrice_Morse& AVV)
745{
746 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
747 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
748 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
749 const IntTab& face_voisins=domaine_VEF.face_voisins();
750 const DoubleTab& normales = domaine_VEF.face_normales();
751 int nb_aretes_tot=domaine_VEF.domaine().nb_aretes();
752 int i, j, k, l;
753 double psc;
754 // Skip vertices -1 which
755 // are at the end of the sommets array:
756 //while (sommets(dplusdeux-1)==-1)
757 // dplusdeux--;
758
759
760
761 int jmax=5;
762 if(elem2==-1) jmax=4;
763 for(i=0; i<3; i++)
764 {
765 int si=sommets[i];
766 for(j=i+1; j<jmax; j++)
767 {
768 int sj=sommets[j];
769 int arete1;
770 if(j<4)
771 arete1 = chercher_arete(domaine_VEF,elem1, si, sj,
772 elem_aretes, aretes_som);
773 else
774 arete1 = chercher_arete(domaine_VEF,elem2, si, sj,
775 elem_aretes, aretes_som);
776 assert(arete1!=-1);
777 if(ok_arete[arete1])
778 {
779 calculer_grad_arete(face, face_voisins, i, j,
780 elem1, elem2,
781 faces_op1[i], faces_op2[i],
782 faces_op1[j], faces_op2[j],
783 normales, gradi);
784 if(arete1<nb_aretes_tot)
785 ARR(arete1,arete1)+=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradi,
786 inverse_quantitee_entrelacee,face);
787 int jj=j;
788 for(k=i; k<3; k++)
789 {
790 int sk=sommets[k];
791 for(l=jj+1; l<jmax; l++)
792 {
793 int sl=sommets[l];
794 int arete2;
795 if(l<4)
796 arete2 = chercher_arete(domaine_VEF,elem1, sl, sk,
797 elem_aretes,
798 aretes_som);
799 else
800 arete2 = chercher_arete(domaine_VEF,elem2, sl, sk,
801 elem_aretes,
802 aretes_som);
803 assert(arete2!=-1);
804 if(ok_arete[arete2])
805 {
806 calculer_grad_arete(face, face_voisins, k, l,
807 elem1, elem2,
808 faces_op1[k], faces_op2[k],
809 faces_op1[l], faces_op2[l],
810 normales, gradj);
811 psc=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradj,
812 inverse_quantitee_entrelacee,face);
813 int tmp=arete1;
814 if(arete1>arete2) swap(arete1, arete2);
815 if(arete1<nb_aretes_tot)
816 if(arete2<nb_aretes_tot)
817 ARR(arete1,arete2)+=psc;
818 else
819 ARV(arete1,arete2-nb_aretes_tot)+=psc;
820 else if(arete2<nb_aretes_tot)
821 AVR(arete1-nb_aretes_tot,arete2)+=psc;
822 else
823 AVV(arete1-nb_aretes_tot,arete2-nb_aretes_tot)+=psc;
824 arete1=tmp;
825 }
826 }
827 jj=k+1;
828 }
829 }
830 }
831 }
832}
833static void
834contribuer_matrice_NeumannP0P1(int elem, const ArrOfInt& sommets, IntLists& voisins, DoubleLists& coeffs)
835
836{
837 int dimension=Objet_U::dimension,
838 dplusun=dimension+1;
839 for(int i=0; i<dplusun; i++)
840 {
841 int si=sommets[i];
842 int rang1=voisins[elem].rang(si);
843 if(rang1==-1)
844 {
845 voisins[elem].add(si);
846 coeffs[elem].add(0);
847 }
848 }
849}
850
851static void
852update_matrice_NeumannP0P1(const Domaine_VEF& domaine_VEF,
853 const DoubleTab& inverse_quantitee_entrelacee,
854 int face, int elem,
855 ArrOfInt& sommets, ArrOfInt& faces_op1, const ArrOfDouble& coef_som,
856 Matrice_Morse& ARR, Matrice_Morse& ARV,
857 Matrice_Morse& AVR, Matrice_Morse& AVV)
858{
859 const DoubleTab& normales = domaine_VEF.face_normales();
860 const IntTab& face_voisins = domaine_VEF.face_voisins();
861
862 int nb_elem_tot=domaine_VEF.nb_elem();
863 int nb_som_tot=domaine_VEF.nb_som();
864
865 int dimension=Objet_U::dimension,
866 dplusun=dimension+1;
867 double unsurdim=1./Objet_U::dimension;
868 double psc;
869 // double coeff_som=1./(dimension)/(dplusun);
870
871
872 for(int i=0; i<dplusun; i++)
873 {
874 int si=sommets[i];
875 calculer_grad(face_voisins, elem, -1, coef_som,si, faces_op1[i],
876 -1, normales, gradi);
877 if(faces_op1[i]!=face)
878 for (int comp=0; comp<dimension; comp++)
879 gradi[comp]+= normales(face,comp)*unsurdim;
880 for(int k=0; k<dimension; k++)
881 gradj[k]=normales(face,k);
882 psc=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradj,
883 inverse_quantitee_entrelacee,face);
884 if(elem<nb_elem_tot)
885 if(si<nb_som_tot)
886 ARR(elem,si)+=psc;
887 else
888 ARV(elem,si-nb_som_tot)+=psc;
889 else if(si<nb_som_tot)
890 AVR(elem-nb_elem_tot,si)+=psc;
891 else
892 AVV(elem-nb_elem_tot,si-nb_som_tot)+=psc;
893 }
894}
895
896static void
897contribuer_matrice_NeumannP1P1(int elem, const ArrOfInt& sommets, IntLists& voisins, DoubleLists& coeffs)
898
899{
900 int dimension=Objet_U::dimension,
901 dplusun=dimension+1;
902 for(int i=0; i<dplusun; i++)
903 {
904 int si=sommets[i];
905 for(int j=i+1; j<dplusun; j++)
906 {
907 int sj=sommets[j];
908 int rang=voisins[si].rang(sj);
909 if(rang==-1)
910 {
911 voisins[si].add(sj);
912 coeffs[si].add(0);
913 }
914 }
915 }
916}
917
918
919static void
920update_matrice_NeumannP1P1(const Domaine_VEF& domaine_VEF,
921 const DoubleTab& inverse_quantitee_entrelacee,
922 int face, int elem,
923 ArrOfInt& sommets, ArrOfInt& faces_op1, const ArrOfDouble& coef_som,
924 Matrice_Morse& ARR, Matrice_Morse& ARV,
925 Matrice_Morse& AVR, Matrice_Morse& AVV)
926{
927 const DoubleTab& normales = domaine_VEF.face_normales();
928 const IntTab& face_voisins = domaine_VEF.face_voisins();
929
930 int dimension=Objet_U::dimension,
931 dplusun=dimension+1;
932 double unsurdim=1./Objet_U::dimension;
933 double psc;
934 // double coeff_som=1./(dimension)/(dplusun);
935 //coeff_som*=coeff_som;
936
937
938 int nb_som_tot=domaine_VEF.nb_som();
939 int i,j;
940 for(i=0; i<dplusun; i++)
941 {
942 int si=sommets[i];
943 calculer_grad(face_voisins, elem, -1, coef_som,si, faces_op1[i],
944 -1, normales, gradi);
945 if(faces_op1[i]!=face)
946 for (int comp=0; comp<dimension; comp++)
947 gradi[comp]+= normales(face,comp)*unsurdim;
948 if(si<nb_som_tot)
949 ARR(si,si)+=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradi,
950 inverse_quantitee_entrelacee,face);
951 for(j=i+1; j<dplusun; j++)
952 {
953 int sj=sommets[j];
954 calculer_grad(face_voisins, elem, -1, coef_som,sj, faces_op1[j],
955 -1, normales, gradj);
956 if(faces_op1[j]!=face)
957 for (int comp=0; comp<dimension; comp++)
958 gradj[comp]+= normales(face,comp)*unsurdim;
959 psc=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradj,inverse_quantitee_entrelacee,face);
960 if(si<nb_som_tot)
961 if(sj<nb_som_tot)
962 ARR(si,sj)+=psc;
963 else
964 ARV(si,sj-nb_som_tot)+=psc;
965 else if(sj<nb_som_tot)
966 AVR(si-nb_som_tot,sj)+=psc;
967 else
968 AVV(si-nb_som_tot,sj-nb_som_tot)+=psc;
969 }
970 }
971}
972
973static void
974contribuer_matrice_SymetrieP1P1(int elem, const ArrOfInt& sommets, IntLists& voisins, DoubleLists& coeffs)
975{
976 int dimension=Objet_U::dimension,
977 dplusun=dimension+1;
978 for(int i=0; i<dplusun; i++)
979 {
980 int si=sommets[i];
981 for(int j=i+1; j<dplusun; j++)
982 {
983 int sj=sommets[j];
984 int rang=voisins[si].rang(sj);
985 if(rang==-1)
986 {
987 voisins[si].add(sj);
988 coeffs[si].add(0);
989 }
990 }
991 }
992}
993
994
995static void
996update_matrice_SymetrieP1P1(const Domaine_VEF& domaine_VEF,
997 const DoubleTab& inverse_quantitee_entrelacee,
998 int face, int elem,
999 ArrOfInt& sommets, ArrOfInt& faces_op1, const ArrOfDouble& coef_som,
1000 Matrice_Morse& ARR, Matrice_Morse& ARV,
1001 Matrice_Morse& AVR, Matrice_Morse& AVV)
1002{
1003 const DoubleTab& normales = domaine_VEF.face_normales();
1004 const IntTab& face_voisins = domaine_VEF.face_voisins();
1005
1006 int dimension=Objet_U::dimension,
1007 dplusun=dimension+1;
1008 double psc;
1009 // double coeff_som=1./(dimension)/(dplusun);
1010 //coeff_som*=coeff_som;
1011
1012
1013 int nb_som_tot=domaine_VEF.nb_som();
1014 int i,j;
1015 for(i=0; i<dplusun; i++)
1016 {
1017 int si=sommets[i];
1018 calculer_grad(face_voisins, elem, -1, coef_som, si, faces_op1[i],
1019 -1, normales, gradi);
1020 projette(gradi, face, normales);
1021 if(si<nb_som_tot)
1022 ARR(si,si)+=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradi,
1023 inverse_quantitee_entrelacee,face);
1024 for(j=i+1; j<dplusun; j++)
1025 {
1026 int sj=sommets[j];
1027 calculer_grad(face_voisins, elem, -1, coef_som,sj, faces_op1[j],
1028 -1, normales, gradj);
1029 projette(gradj, face, normales);
1030 psc=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradj,
1031 inverse_quantitee_entrelacee,face);
1032 if(si<nb_som_tot)
1033 if(sj<nb_som_tot)
1034 ARR(si,sj)+=psc;
1035 else
1036 ARV(si,sj-nb_som_tot)+=psc;
1037 else if(sj<nb_som_tot)
1038 AVR(si-nb_som_tot,sj)+=psc;
1039 else
1040 AVV(si-nb_som_tot,sj-nb_som_tot)+=psc;
1041 }
1042 }
1043}
1044
1045static void
1046contribuer_matrice_NeumannPaPa(const Domaine_VEF& domaine_VEF, int elem, const ArrOfInt& sommets, IntLists& voisins, DoubleLists& coeffs)
1047
1048{
1049 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1050 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1051 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1052 for(int i=0; i<3; i++)
1053 {
1054 int si=sommets[i];
1055 for(int j=i+1; j<4; j++)
1056 {
1057 int sj=sommets[j];
1058 int arete1;
1059 arete1 = chercher_arete(domaine_VEF,elem, si, sj,
1060 elem_aretes, aretes_som);
1061 if(ok_arete[arete1])
1062 {
1063 int jj=j;
1064 for(int k=i; k<3; k++)
1065 {
1066 int sk=sommets[k];
1067 for(int l=jj+1; l<4; l++)
1068 {
1069 int sl=sommets[l];
1070 int arete2= chercher_arete(domaine_VEF,elem, sl, sk,
1071 elem_aretes,
1072 aretes_som);
1073 assert(arete2!=-1);
1074 if(ok_arete[arete2])
1075 {
1076 int tmp=arete1;
1077 if(arete1>arete2) swap(arete1, arete2);
1078 int rang=voisins[arete1].rang(arete2);
1079 if(rang==-1)
1080 {
1081 voisins[arete1].add(arete2);
1082 coeffs[arete1].add(0);
1083
1084 }
1085 arete1=tmp;
1086 }
1087 }
1088 jj=k+1;
1089 }
1090 }
1091 }
1092 }
1093}
1094
1095static void
1096update_matrice_NeumannPaPa(const Domaine_VEF& domaine_VEF,
1097 const DoubleTab& inverse_quantitee_entrelacee,
1098 int face, int elem,
1099 ArrOfInt& sommets, ArrOfInt& faces_op1,
1100 Matrice_Morse& ARR, Matrice_Morse& ARV,
1101 Matrice_Morse& AVR, Matrice_Morse& AVV)
1102{
1103 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1104 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1105 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1106 const IntTab& face_voisins=domaine_VEF.face_voisins();
1107 const DoubleTab& normales = domaine_VEF.face_normales();
1108 int nb_aretes_tot=domaine_VEF.domaine().nb_aretes();
1109
1110 int i, j, k, l;
1111 double psc;
1112
1113
1114 for(i=0; i<3; i++)
1115 {
1116 int si=sommets[i];
1117 for(j=i+1; j<4; j++)
1118 {
1119 int sj=sommets[j];
1120 int arete1;
1121 arete1 = chercher_arete(domaine_VEF,elem, si, sj,
1122 elem_aretes, aretes_som);
1123 if(ok_arete[arete1])
1124 {
1125 calculer_grad_arete(face, face_voisins, i, j,
1126 elem, -1,
1127 faces_op1[i], -1,
1128 faces_op1[j], -1,
1129 normales, gradi);
1130 if(arete1<nb_aretes_tot)
1131 ARR(arete1,arete1)+=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradi,
1132 inverse_quantitee_entrelacee,face);
1133 int jj=j;
1134 for(k=i; k<3; k++)
1135 {
1136 int sk=sommets[k];
1137 for(l=jj+1; l<4; l++)
1138 {
1139 int sl=sommets[l];
1140 int arete2= chercher_arete(domaine_VEF,elem, sl, sk,
1141 elem_aretes,
1142 aretes_som);
1143 assert(arete2!=-1);
1144 if(ok_arete[arete2])
1145 {
1146 calculer_grad_arete(face, face_voisins, k, l,
1147 elem, -1,
1148 faces_op1[k], -1,
1149 faces_op1[l], -1,
1150 normales, gradj);
1151 psc=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradj,
1152 inverse_quantitee_entrelacee,face);
1153 int tmp=arete1;
1154 if(arete1>arete2) swap(arete1, arete2);
1155 if(arete1<nb_aretes_tot)
1156 if(arete2<nb_aretes_tot)
1157 ARR(arete1,arete2)+=psc;
1158 else
1159 ARV(arete1,arete2-nb_aretes_tot)+=psc;
1160 else if(arete2<nb_aretes_tot)
1161 AVR(arete1-nb_aretes_tot,arete2)+=psc;
1162 else
1163 AVV(arete1-nb_aretes_tot,arete2-nb_aretes_tot)+=psc;
1164 arete1=tmp;
1165 }
1166 }
1167 jj=k+1;
1168 }
1169 }
1170 }
1171 }
1172}
1173
1174static void
1175contribuer_matrice_SymetriePaPa(const Domaine_VEF& domaine_VEF, int elem, const ArrOfInt& sommets, IntLists& voisins, DoubleLists& coeffs)
1176{
1177 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1178 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1179 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1180 for(int i=0; i<3; i++)
1181 {
1182 int si=sommets[i];
1183 for(int j=i+1; j<4; j++)
1184 {
1185 int sj=sommets[j];
1186 int arete1;
1187 arete1 = chercher_arete(domaine_VEF,elem, si, sj,
1188 elem_aretes, aretes_som);
1189 if(ok_arete[arete1])
1190 {
1191 int jj=j;
1192 for(int k=i; k<3; k++)
1193 {
1194 int sk=sommets[k];
1195 for(int l=jj+1; l<4; l++)
1196 {
1197 int sl=sommets[l];
1198 int arete2= chercher_arete(domaine_VEF,elem, sl, sk,
1199 elem_aretes,
1200 aretes_som);
1201 assert(arete2!=-1);
1202 if(ok_arete[arete2])
1203 {
1204 int tmp=arete1;
1205 if(arete1>arete2) swap(arete1, arete2);
1206 int rang=voisins[arete1].rang(arete2);
1207 if(rang==-1)
1208 {
1209 voisins[arete1].add(arete2);
1210 coeffs[arete1].add(0);
1211
1212 }
1213 arete1=tmp;
1214 }
1215 }
1216 jj=k+1;
1217 }
1218 }
1219 }
1220 }
1221}
1222
1223static void
1224update_matrice_SymetriePaPa(const Domaine_VEF& domaine_VEF,
1225 const DoubleTab& inverse_quantitee_entrelacee,
1226 int face, int elem,
1227 ArrOfInt& sommets, ArrOfInt& faces_op1,
1228 Matrice_Morse& ARR, Matrice_Morse& ARV,
1229 Matrice_Morse& AVR, Matrice_Morse& AVV)
1230{
1231 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1232 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1233 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1234 const IntTab& face_voisins=domaine_VEF.face_voisins();
1235 const DoubleTab& normales = domaine_VEF.face_normales();
1236 int i, j, k, l;
1237 int nb_aretes_tot=domaine_VEF.domaine().nb_aretes();
1238
1239 double psc;
1240
1241
1242 for(i=0; i<3; i++)
1243 {
1244 int si=sommets[i];
1245 for(j=i+1; j<4; j++)
1246 {
1247 int sj=sommets[j];
1248 int arete1;
1249 arete1 = chercher_arete(domaine_VEF,elem, si, sj,
1250 elem_aretes, aretes_som);
1251 if(ok_arete[arete1])
1252 {
1253 calculer_grad_arete(face, face_voisins, i, j,
1254 elem, -1,
1255 faces_op1[i], -1,
1256 faces_op1[j], -1,
1257 normales, gradi);
1258 projette(gradi, face, normales);
1259 if(arete1<nb_aretes_tot)
1260 ARR(arete1,arete1)+=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradi,
1261 inverse_quantitee_entrelacee,face);
1262 int jj=j;
1263 for(k=i; k<3; k++)
1264 {
1265 int sk=sommets[k];
1266 for(l=jj+1; l<4; l++)
1267 {
1268 int sl=sommets[l];
1269 int arete2= chercher_arete(domaine_VEF,elem, sl, sk,
1270 elem_aretes,
1271 aretes_som);
1272 assert(arete2!=-1);
1273 if(ok_arete[arete2])
1274 {
1275 calculer_grad_arete(face, face_voisins, k, l,
1276 elem, -1,
1277 faces_op1[k], -1,
1278 faces_op1[l], -1,
1279 normales, gradj);
1280 projette(gradj, face, normales);
1281 psc=dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradj,
1282 inverse_quantitee_entrelacee,face);
1283 int tmp=arete1;
1284 if(arete1>arete2) swap(arete1, arete2);
1285 if(arete1<nb_aretes_tot)
1286 if(arete2<nb_aretes_tot)
1287 ARR(arete1,arete2)+=psc;
1288 else
1289 ARV(arete1,arete2-nb_aretes_tot)+=psc;
1290 else if(arete2<nb_aretes_tot)
1291 AVR(arete1-nb_aretes_tot,arete2)+=psc;
1292 else
1293 AVV(arete1-nb_aretes_tot,arete2-nb_aretes_tot)+=psc;
1294 arete1=tmp;
1295 }
1296 }
1297 jj=k+1;
1298 }
1299 }
1300 }
1301 }
1302}
1303
1304static void contribuer_matriceP0Pa(const Domaine_VEF& domaine_VEF, int elem1, int elem2,
1305 const ArrOfInt& sommets,
1306 IntLists& voisins,
1307 DoubleLists& coeffs)
1308
1309{
1310 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1311 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1312 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1313 int jmax=5;
1314 if(elem2==-1) jmax=4;
1315 for(int i=0; i<3; i++)
1316 {
1317 int si=sommets[i];
1318 for(int j=i+1; j<jmax; j++)
1319 {
1320 int sj=sommets[j];
1321 int arete1;
1322 if(j<4)
1323 arete1 = chercher_arete(domaine_VEF,elem1, si, sj,
1324 elem_aretes, aretes_som);
1325 else
1326 arete1 = chercher_arete(domaine_VEF,elem2, si, sj,
1327 elem_aretes, aretes_som);
1328 if(ok_arete[arete1])
1329 {
1330 int rang=voisins[elem1].rang(arete1);
1331 if(rang==-1)
1332 {
1333 voisins[elem1].add(arete1);
1334 coeffs[elem1].add(0);
1335 }
1336 if(elem2!=-1)
1337 {
1338 int rangbis=voisins[elem2].rang(arete1);
1339 if(rangbis==-1)
1340 {
1341 voisins[elem2].add(arete1);
1342 coeffs[elem2].add(0);
1343 }
1344 }
1345 }
1346 }
1347 }
1348}
1349
1350static void update_matriceP0Pa(const Domaine_VEF& domaine_VEF,
1351 const DoubleTab& inverse_quantitee_entrelacee,
1352 int face, int elem1, int elem2,
1353 ArrOfInt& sommets,
1354 ArrOfInt& faces_op1,
1355 ArrOfInt& faces_op2,
1356 Matrice_Morse& ARR, Matrice_Morse& ARV,
1357 Matrice_Morse& AVR, Matrice_Morse& AVV)
1358{
1359 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1360 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1361 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1362 const IntTab& face_voisins=domaine_VEF.face_voisins();
1363 const DoubleTab& normales = domaine_VEF.face_normales();
1364 int nb_aretes_tot=domaine_VEF.domaine().nb_aretes();
1365 int nb_elem_tot=domaine_VEF.domaine().nb_elem();
1366
1367 int i, j;
1368 double psc;
1369
1370
1371 int jmax=5;
1372 if(elem2==-1) jmax=4;
1373 for(i=0; i<3; i++)
1374 {
1375 int si=sommets[i];
1376 for(j=i+1; j<jmax; j++)
1377 {
1378 int sj=sommets[j];
1379 int arete1;
1380 if(j<4)
1381 arete1 = chercher_arete(domaine_VEF,elem1, si, sj,
1382 elem_aretes, aretes_som);
1383 else
1384 arete1 = chercher_arete(domaine_VEF,elem2, si, sj,
1385 elem_aretes, aretes_som);
1386 if(ok_arete[arete1])
1387 {
1388 calculer_grad_arete(face, face_voisins, i, j,
1389 elem1, elem2,
1390 faces_op1[i], faces_op2[i],
1391 faces_op1[j], faces_op2[j],
1392 normales, gradi);
1393 psc=0;
1394 for(int comp=0; comp<3; comp++)
1395 psc+=gradi[comp]*normales(face, comp)
1396 *(-inverse_quantitee_entrelacee(face,comp));
1397 if(elem1<nb_elem_tot)
1398 if(arete1<nb_aretes_tot)
1399 ARR(elem1,arete1)+=psc;
1400 else
1401 ARV(elem1,arete1-nb_aretes_tot)+=psc;
1402 else if(arete1<nb_aretes_tot)
1403 AVR(elem1-nb_elem_tot,arete1)+=psc;
1404 else
1405 AVV(elem1-nb_elem_tot,arete1-nb_aretes_tot)+=psc;
1406
1407 if(elem2!=-1)
1408 {
1409 psc*=-1.0;
1410 if(elem2<nb_elem_tot)
1411 if(arete1<nb_aretes_tot)
1412 ARR(elem2,arete1)+=psc;
1413 else
1414 ARV(elem2,arete1-nb_aretes_tot)+=psc;
1415 else if(arete1<nb_aretes_tot)
1416 AVR(elem2-nb_elem_tot,arete1)+=psc;
1417 else
1418 AVV(elem2-nb_elem_tot,arete1-nb_aretes_tot)+=psc;
1419 }
1420 }
1421 }
1422 }
1423}
1424
1425static void
1426contribuer_matrice_NeumannP0Pa(const Domaine_VEF& domaine_VEF, int elem,
1427 const ArrOfInt& sommets,
1428 IntLists& voisins,
1429 DoubleLists& coeffs)
1430
1431{
1432 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1433 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1434 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1435 for(int i=0; i<3; i++)
1436 {
1437 int si=sommets[i];
1438 for(int j=i+1; j<4; j++)
1439 {
1440 int sj=sommets[j];
1441 int arete1;
1442 arete1 = chercher_arete(domaine_VEF,elem, si, sj,
1443 elem_aretes, aretes_som);
1444 if(ok_arete[arete1])
1445 {
1446 int rang=voisins[elem].rang(arete1);
1447 if(rang==-1)
1448 {
1449 voisins[elem].add(arete1);
1450 coeffs[elem].add(0);
1451 }
1452 }
1453 }
1454 }
1455}
1456
1457static void
1458update_matrice_NeumannP0Pa(const Domaine_VEF& domaine_VEF,
1459 const DoubleTab& inverse_quantitee_entrelacee,
1460 int face, int elem,
1461 ArrOfInt& sommets,
1462 ArrOfInt& faces_op1,
1463 Matrice_Morse& ARR, Matrice_Morse& ARV,
1464 Matrice_Morse& AVR, Matrice_Morse& AVV)
1465{
1466 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1467 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1468 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1469 const IntTab& face_voisins=domaine_VEF.face_voisins();
1470 const DoubleTab& normales = domaine_VEF.face_normales();
1471 int nb_aretes_tot=domaine_VEF.domaine().nb_aretes();
1472 int nb_elem_tot=domaine_VEF.domaine().nb_elem();
1473
1474 int i, j;
1475 double psc;
1476
1477
1478 for(i=0; i<3; i++)
1479 {
1480 int si=sommets[i];
1481 for(j=i+1; j<4; j++)
1482 {
1483 int sj=sommets[j];
1484 int arete1;
1485 arete1 = chercher_arete(domaine_VEF,elem, si, sj,
1486 elem_aretes, aretes_som);
1487 if(ok_arete[arete1])
1488 {
1489 calculer_grad_arete(face, face_voisins, i, j,
1490 elem, -1,
1491 faces_op1[i], -1,
1492 faces_op1[j], -1,
1493 normales, gradi);
1494 psc=0;
1495 for(int comp=0; comp<3; comp++)
1496 psc+=gradi[comp]*normales(face, comp)
1497 *(-inverse_quantitee_entrelacee(face,comp));
1498 if(elem<nb_elem_tot)
1499 if(arete1<nb_aretes_tot)
1500 ARR(elem,arete1)+=psc;
1501 else
1502 ARV(elem,arete1-nb_aretes_tot)+=psc;
1503 else if(arete1<nb_aretes_tot)
1504 AVR(elem-nb_elem_tot,arete1)+=psc;
1505 else
1506 AVV(elem-nb_elem_tot,arete1-nb_aretes_tot)+=psc;
1507 }
1508 }
1509 }
1510}
1511
1512static void contribuer_matriceP1Pa(const Domaine_VEF& domaine_VEF, int elem1, int elem2,
1513 const ArrOfInt& sommets,
1514 IntLists& voisins,
1515 DoubleLists& coeffs)
1516
1517{
1518 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1519 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1520 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1521 int jmax=5;
1522 if(elem2==-1) jmax=4;
1523 for(int i=0; i<3; i++)
1524 {
1525 int si=sommets[i];
1526 for(int j=i+1; j<jmax; j++)
1527 {
1528 int sj=sommets[j];
1529 int arete1;
1530 if(j<4)
1531 arete1 = chercher_arete(domaine_VEF,elem1, si, sj,
1532 elem_aretes, aretes_som);
1533 else
1534 arete1 = chercher_arete(domaine_VEF,elem2, si, sj,
1535 elem_aretes, aretes_som);
1536 if(ok_arete[arete1])
1537 {
1538 for(int k=0; k<jmax; k++)
1539 {
1540 int sk=sommets[k];
1541 int rang1=voisins[sk].rang(arete1);
1542 if(rang1==-1)
1543 {
1544 voisins[sk].add(arete1);
1545 coeffs[sk].add(0);
1546 }
1547 }
1548 }
1549 }
1550 }
1551}
1552
1553static void update_matriceP1Pa(const Domaine_VEF& domaine_VEF,
1554 const DoubleTab& inverse_quantitee_entrelacee,
1555 int face, int elem1, int elem2,
1556 ArrOfInt& sommets,
1557 ArrOfInt& faces_op1,
1558 ArrOfInt& faces_op2, const ArrOfDouble& coef_som,
1559 Matrice_Morse& ARR, Matrice_Morse& ARV,
1560 Matrice_Morse& AVR, Matrice_Morse& AVV)
1561{
1562 // int dimension=Objet_U::dimension,
1563 // dplusun=dimension+1;
1564 // double coeff_som=1./(dimension)/(dplusun);
1565 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1566 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1567 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1568 const IntTab& face_voisins=domaine_VEF.face_voisins();
1569 const DoubleTab& normales = domaine_VEF.face_normales();
1570 int nb_aretes_tot=domaine_VEF.domaine().nb_aretes();
1571 int nb_som_tot=domaine_VEF.domaine().nb_som();
1572
1573 int i, j, k;
1574 double psc;
1575
1576
1577 int jmax=5;
1578 if(elem2==-1) jmax=4;
1579 for(i=0; i<3; i++)
1580 {
1581 int si=sommets[i];
1582 for(j=i+1; j<jmax; j++)
1583 {
1584 int sj=sommets[j];
1585 int arete1;
1586 if(j<4)
1587 arete1 = chercher_arete(domaine_VEF,elem1, si, sj,
1588 elem_aretes, aretes_som);
1589 else
1590 arete1 = chercher_arete(domaine_VEF,elem2, si, sj,
1591 elem_aretes, aretes_som);
1592 if(ok_arete[arete1])
1593 {
1594 calculer_grad_arete(face, face_voisins, i, j,
1595 elem1, elem2,
1596 faces_op1[i], faces_op2[i],
1597 faces_op1[j], faces_op2[j],
1598 normales, gradi);
1599 for(k=0; k<jmax; k++)
1600 {
1601 int sk=sommets[k];
1602 calculer_grad(face_voisins, elem1, elem2, coef_som, sk,
1603 faces_op1[k], faces_op2[k],
1604 normales, gradj);
1605 psc=-dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradj,
1606 inverse_quantitee_entrelacee,face);
1607 if(sk<nb_som_tot)
1608 if(arete1<nb_aretes_tot)
1609 ARR(sk,arete1)+=psc;
1610 else
1611 ARV(sk,arete1-nb_aretes_tot)+=psc;
1612 else if(arete1<nb_aretes_tot)
1613 AVR(sk-nb_som_tot,arete1)+=psc;
1614 else
1615 AVV(sk-nb_som_tot,arete1-nb_aretes_tot)+=psc;
1616 }
1617 }
1618 }
1619 }
1620}
1621
1622static void
1623contribuer_matrice_NeumannP1Pa(const Domaine_VEF& domaine_VEF, int elem,
1624 const ArrOfInt& sommets,
1625 IntLists& voisins,
1626 DoubleLists& coeffs)
1627
1628{
1629 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1630 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1631 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1632 for(int i=0; i<3; i++)
1633 {
1634 int si=sommets[i];
1635 for(int j=i+1; j<4; j++)
1636 {
1637 int sj=sommets[j];
1638 int arete1;
1639 arete1 = chercher_arete(domaine_VEF,elem, si, sj,
1640 elem_aretes, aretes_som);
1641 if(ok_arete[arete1])
1642 {
1643 for(int k=0; k<4; k++)
1644 {
1645 int sk=sommets[k];
1646 int rang1=voisins[sk].rang(arete1);
1647 if(rang1==-1)
1648 {
1649 voisins[sk].add(arete1);
1650 coeffs[sk].add(0);
1651
1652 }
1653 }
1654 }
1655 }
1656 }
1657}
1658
1659static void
1660update_matrice_NeumannP1Pa(const Domaine_VEF& domaine_VEF,
1661 const DoubleTab& inverse_quantitee_entrelacee,
1662 int face, int elem,
1663 ArrOfInt& sommets,
1664 ArrOfInt& faces_op1, const ArrOfDouble& coef_som,
1665 Matrice_Morse& ARR, Matrice_Morse& ARV,
1666 Matrice_Morse& AVR, Matrice_Morse& AVV)
1667{
1668 int dimension=Objet_U::dimension;
1669 // dplusun=dimension+1;
1670 // double coeff_som=1./(dimension)/(dplusun);
1671 double unsurdim=1./Objet_U::dimension;
1672 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1673 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1674 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1675 const IntTab& face_voisins=domaine_VEF.face_voisins();
1676 const DoubleTab& normales = domaine_VEF.face_normales();
1677 int nb_aretes_tot=domaine_VEF.domaine().nb_aretes();
1678 int nb_som_tot=domaine_VEF.domaine().nb_som();
1679
1680 int i, j, k;
1681 double psc;
1682
1683
1684 for(i=0; i<3; i++)
1685 {
1686 int si=sommets[i];
1687 for(j=i+1; j<4; j++)
1688 {
1689 int sj=sommets[j];
1690 int arete1;
1691 arete1 = chercher_arete(domaine_VEF,elem, si, sj,
1692 elem_aretes, aretes_som);
1693 if(ok_arete[arete1])
1694 {
1695 calculer_grad_arete(face, face_voisins, i, j,
1696 elem, -1,
1697 faces_op1[i], -1,
1698 faces_op1[j], -1,
1699 normales, gradi);
1700 for(k=0; k<4; k++)
1701 {
1702 int sk=sommets[k];
1703 calculer_grad(face_voisins, elem, -1, coef_som,sk,
1704 faces_op1[k], -1,
1705 normales, gradj);
1706 if(faces_op1[k]!=face)
1707 for (int comp=0; comp<dimension; comp++)
1708 gradj[comp]+= normales(face,comp)*unsurdim;
1709 psc=-dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradj,
1710 inverse_quantitee_entrelacee,face);
1711 if(sk<nb_som_tot)
1712 if(arete1<nb_aretes_tot)
1713 ARR(sk,arete1)+=psc;
1714 else
1715 ARV(sk,arete1-nb_aretes_tot)+=psc;
1716 else if(arete1<nb_aretes_tot)
1717 AVR(sk-nb_som_tot,arete1)+=psc;
1718 else
1719 AVV(sk-nb_som_tot,arete1-nb_aretes_tot)+=psc;
1720 }
1721 }
1722 }
1723 }
1724}
1725
1726static void
1727contribuer_matrice_SymetrieP1Pa(const Domaine_VEF& domaine_VEF, int elem,
1728 const ArrOfInt& sommets,
1729 IntLists& voisins,
1730 DoubleLists& coeffs)
1731
1732{
1733 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1734 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1735 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1736 for(int i=0; i<3; i++)
1737 {
1738 int si=sommets[i];
1739 for(int j=i+1; j<4; j++)
1740 {
1741 int sj=sommets[j];
1742 int arete1;
1743 arete1 = chercher_arete(domaine_VEF,elem, si, sj,
1744 elem_aretes, aretes_som);
1745 if(ok_arete[arete1])
1746 {
1747 for(int k=0; k<4; k++)
1748 {
1749 int sk=sommets[k];
1750 int rang1=voisins[sk].rang(arete1);
1751 if(rang1==-1)
1752 {
1753 voisins[sk].add(arete1);
1754 coeffs[sk].add(0);
1755 }
1756 }
1757 }
1758 }
1759 }
1760}
1761
1762static void
1763update_matrice_SymetrieP1Pa(const Domaine_VEF& domaine_VEF,
1764 const DoubleTab& inverse_quantitee_entrelacee,
1765 int face, int elem,
1766 ArrOfInt& sommets,
1767 ArrOfInt& faces_op1, const ArrOfDouble& coef_som,
1768 Matrice_Morse& ARR, Matrice_Morse& ARV,
1769 Matrice_Morse& AVR, Matrice_Morse& AVV)
1770{
1771 // int dimension=Objet_U::dimension,
1772 // dplusun=dimension+1;
1773 // double coeff_som=1./(dimension)/(dplusun);
1774 const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
1775 const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
1776 const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
1777 const IntTab& face_voisins=domaine_VEF.face_voisins();
1778 const DoubleTab& normales = domaine_VEF.face_normales();
1779 int nb_aretes_tot=domaine_VEF.domaine().nb_aretes();
1780 int nb_som_tot=domaine_VEF.domaine().nb_som();
1781
1782 int i, j, k;
1783 double psc;
1784
1785
1786 for(i=0; i<3; i++)
1787 {
1788 int si=sommets[i];
1789 for(j=i+1; j<4; j++)
1790 {
1791 int sj=sommets[j];
1792 int arete1;
1793 arete1 = chercher_arete(domaine_VEF,elem, si, sj,
1794 elem_aretes, aretes_som);
1795 if(ok_arete[arete1])
1796 {
1797 calculer_grad_arete(face, face_voisins, i, j,
1798 elem, -1,
1799 faces_op1[i], -1,
1800 faces_op1[j], -1,
1801 normales, gradi);
1802 projette(gradi, face, normales);
1803 for(k=0; k<4; k++)
1804 {
1805 int sk=sommets[k];
1806 calculer_grad(face_voisins, elem, -1, coef_som,sk,
1807 faces_op1[k], -1,
1808 normales, gradj);
1809 projette(gradj, face, normales);
1810 psc=-dotproduct_array_fois_inverse_quantitee_entrelacee(gradi,gradj,
1811 inverse_quantitee_entrelacee,face);
1812 if(sk<nb_som_tot)
1813 if(arete1<nb_aretes_tot)
1814 ARR(sk,arete1)+=psc;
1815 else
1816 ARV(sk,arete1-nb_aretes_tot)+=psc;
1817 else if(arete1<nb_aretes_tot)
1818 AVR(sk-nb_som_tot,arete1)+=psc;
1819 else
1820 AVV(sk-nb_som_tot,arete1-nb_aretes_tot)+=psc;
1821 }
1822 }
1823 }
1824 }
1825}
1826
1827void
1828assemblerP0P0(const Domaine_dis_base& z,
1829 const Domaine_Cl_dis_base& zcl,
1830 Matrice& matrice,
1831 const DoubleTab& inverse_quantitee_entrelacee)
1832{
1833 Assembleur_P_VEF Assembleur_P0;
1834 Assembleur_P0.associer_domaine_dis_base(z);
1835 Assembleur_P0.associer_domaine_cl_dis_base(zcl);
1836 Assembleur_P0.remplir(matrice,inverse_quantitee_entrelacee);
1837 Cerr << "Assemblage P0 OK" << finl;
1838}
1839
1840void
1841updateP0P0(const Domaine_dis_base& z,
1842 const Domaine_Cl_dis_base& zcl,
1843 Matrice& matrice,
1844 const DoubleTab& inverse_quantitee_entrelacee)
1845{
1846 Assembleur_P_VEF Assembleur_P0;
1847 Assembleur_P0.associer_domaine_dis_base(z);
1848 Assembleur_P0.associer_domaine_cl_dis_base(zcl);
1849 Assembleur_P0.remplir(matrice,inverse_quantitee_entrelacee);
1850 Cerr << "Update P0 OK" << finl;
1851}
1852
1853void assemblerP1P1(const Domaine_dis_base& z,
1854 const Domaine_Cl_dis_base& zcl,
1855 Matrice& matrice, const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som)
1856{
1857 int dimension=Objet_U::dimension;
1858 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, z);
1859 const Domaine& domaine=domaine_VEF.domaine();
1860 const Domaine_Cl_VEF& domaine_Cl_VEF=ref_cast(Domaine_Cl_VEF, zcl);
1861 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
1862 const IntTab& face_voisins = domaine_VEF.face_voisins();
1863 int nint = domaine_VEF.premiere_face_int();
1864 int nb_faces = domaine_VEF.nb_faces_tot();
1865 int nb_som = domaine_VEF.domaine().nb_som_tot();
1866 int elem1, elem2, face, ok;
1867 IntLists voisins(nb_som);
1868 DoubleLists coeffs(nb_som);
1869 ArrOfInt sommets(dimension+2);
1870 ArrOfInt face_opp1(dimension+2);
1871 ArrOfInt face_opp2(dimension+2);
1872 // Boundary faces:
1873 for(int i=0; i<les_cl.size(); i++)
1874 {
1875 const Cond_lim& la_cl = les_cl[i];
1876 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1877 int nb_faces_bord_tot = le_bord.nb_faces_tot();
1878 for(int ind_face=0; ind_face<nb_faces_bord_tot; ind_face++)
1879 {
1880 ok=okface(ind_face, face, la_cl);
1881 if (ok==-1) break;
1882 elem1=face_voisins(face, 0);
1883 elem2=face_voisins(face, 1);
1884 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets,
1885 face_opp1, face_opp2);
1886 sort(sommets, face_opp1, face_opp2);
1887 if(ok==3)
1888 {
1889 contribuer_matrice_NeumannP1P1(elem1, sommets, voisins, coeffs);
1890 }
1891 else if(ok==4)
1892 {
1893 contribuer_matrice_SymetrieP1P1(elem1, sommets, voisins, coeffs);
1894 }
1895 else
1896 contribuer_matriceP1P1(elem1, elem2, sommets, voisins, coeffs);
1897 }
1898 }
1899 face_associee=-1;
1900 for(face=nint; face<nb_faces; face++)
1901 {
1902 elem1=face_voisins(face, 0);
1903 elem2=face_voisins(face, 1);
1904 if (!domaine_VEF.est_une_face_virt_bord(face)) // Process internal faces only
1905 {
1906 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
1907 sort(sommets, face_opp1, face_opp2);
1908 contribuer_matriceP1P1(elem1, elem2, sommets, voisins, coeffs);
1909 }
1910 }
1911 DoubleVect diag(nb_som);
1912 diag=1;
1913 matrice.typer("Matrice_Bloc");
1914 Matrice_Bloc& matrice_bloc=ref_cast(Matrice_Bloc, matrice.valeur());
1915 matrice_bloc.remplir(voisins, coeffs, diag, domaine.nb_som(), domaine.nb_som_tot());
1916 Cerr << "Assemblage P1 OK" << finl;
1917}
1918
1919
1920void updateP1P1(const Domaine_dis_base& z,
1921 const Domaine_Cl_dis_base& zcl,
1922 Matrice& matrice,
1923 const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som)
1924{
1925 int dimension=Objet_U::dimension;
1926 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, z);
1927 const Domaine_Cl_VEF& domaine_Cl_VEF=ref_cast(Domaine_Cl_VEF, zcl);
1928 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
1929 const IntTab& face_voisins = domaine_VEF.face_voisins();
1930 int nint = domaine_VEF.premiere_face_int();
1931 int nb_faces = domaine_VEF.nb_faces_tot();
1932 int elem1, elem2, face, ok;
1933 ArrOfInt sommets(dimension+2);
1934 ArrOfInt face_opp1(dimension+2);
1935 ArrOfInt face_opp2(dimension+2);
1936 Matrice_Bloc& A=ref_cast(Matrice_Bloc, matrice.valeur());
1937 Matrice_Morse_Sym& ARR=ref_cast(Matrice_Morse_Sym, A.get_bloc(0,0).valeur());
1938 Matrice_Morse& ARV=ref_cast(Matrice_Morse, A.get_bloc(0,1).valeur());
1939 Matrice_Morse& AVR=ref_cast(Matrice_Morse, A.get_bloc(1,0).valeur());
1940 Matrice_Morse& AVV=ref_cast(Matrice_Morse, A.get_bloc(1,1).valeur());
1941 // Boundary faces:
1942 for (auto &itr : les_cl)
1943 {
1944 const Cond_lim& la_cl = itr;
1945 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
1946 int nb_faces_bord_tot = le_bord.nb_faces_tot();
1947 for (int ind_face = 0; ind_face < nb_faces_bord_tot; ind_face++)
1948 {
1949 ok = okface(ind_face, face, la_cl);
1950 if (ok == -1)
1951 break;
1952 elem1 = face_voisins(face, 0);
1953 elem2 = face_voisins(face, 1);
1954 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
1955 sort(sommets, face_opp1, face_opp2);
1956 if (ok == 3)
1957 update_matrice_NeumannP1P1(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, sommets, face_opp1, coef_som, ARR, ARV, AVR, AVV);
1958 else if (ok == 4)
1959 update_matrice_SymetrieP1P1(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, sommets, face_opp1, coef_som, ARR, ARV, AVR, AVV);
1960 else
1961 update_matriceP1P1(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, elem2, sommets, face_opp1, face_opp2, coef_som, ARR, ARV, AVR, AVV);
1962 }
1963 }
1964 face_associee=-1;
1965 for (face = nint; face < nb_faces; face++)
1966 {
1967 elem1 = face_voisins(face, 0);
1968 elem2 = face_voisins(face, 1);
1969 if (!domaine_VEF.est_une_face_virt_bord(face)) // Process internal faces only
1970 {
1971 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
1972 sort(sommets, face_opp1, face_opp2);
1973 update_matriceP1P1(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, elem2, sommets, face_opp1, face_opp2, coef_som, ARR, ARV, AVR, AVV);
1974 }
1975 }
1976 int nb_som=domaine_VEF.domaine().nb_som();
1977 for(int i=0; i<nb_som; i++)
1978 if(ARR(i,i)==0)
1979 {
1980 //Cerr << "Modifying orphan vertex row " << i << finl;
1981 ARR(i,i)=1.;
1982 }
1983 Cerr << "Update P1 OK" << finl;
1984}
1985
1986
1987void modifieP1P1neumann(const Domaine_dis_base& z,
1988 const Domaine_Cl_dis_base& zcl,
1989 Matrice& matrice,
1990 const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som)
1991{
1992
1993 //int dimension=Objet_U::dimension;
1994 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, z);
1995 // const Domaine& domaine=domaine_VEF.domaine();
1996 const Domaine_Cl_VEF& domaine_Cl_VEF=ref_cast(Domaine_Cl_VEF, zcl);
1997 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
1998
1999 // int nnz=nb_som;
2000 //ArrOfInt sommets(dimension+2);
2001 //ArrOfInt face_opp1(dimension+2);
2002 //ArrOfInt face_opp2(dimension+2);
2003 Matrice_Bloc& A=ref_cast(Matrice_Bloc, matrice.valeur());
2004 Matrice_Morse_Sym& ARR=ref_cast(Matrice_Morse_Sym, A.get_bloc(0,0).valeur());
2005 // Boundary faces:
2006 assert(ref_cast(Domaine_VEF, z).get_cl_pression_sommet_faible()==0);
2007 int nb_som_tot=z.nb_som();
2008 for(auto& itr : les_cl)
2009 {
2010
2011 const Cond_lim_base& la_cl = itr.valeur();
2012 if (sub_type(Neumann_sortie_libre,la_cl))
2013 {
2014 const Front_VF& le_bord = ref_cast(Front_VF,la_cl.frontiere_dis());
2015 int nb_faces_bord = le_bord.nb_faces_tot();
2016 const IntTab& faces=domaine_VEF.face_sommets();
2017 int nbsf=faces.dimension(1);
2018 for(int ind_face=0; ind_face<nb_faces_bord; ind_face++)
2019 {
2020 int face=le_bord.num_face(ind_face);
2021 for(int som=0; som<nbsf; som++)
2022 {
2023
2024 int som_glob=faces(face,som);
2025 if (som_glob<nb_som_tot)
2026 ARR(som_glob,som_glob)=1e12;
2027 // Cout<<ref_cast(Domaine_VEF, z).numero_premier_sommet()<<" ici "<<som_glob<<finl;
2028 }
2029 }
2030 }
2031 }
2032 Cerr << "Modifie P1P1 Neumann OK" << finl;
2033}
2034
2035void assemblerPaPa(const Domaine_dis_base& z,
2036 const Domaine_Cl_dis_base& zcl,
2037 Matrice& matrice,
2038 const DoubleTab& inverse_quantitee_entrelacee)
2039{
2040 int dimension=Objet_U::dimension;
2041 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, z);
2042 const Domaine& domaine=domaine_VEF.domaine();
2043 const Domaine_Cl_VEF& domaine_Cl_VEF=ref_cast(Domaine_Cl_VEF, zcl);
2044 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
2045 const IntTab& face_voisins = domaine_VEF.face_voisins();
2046 int nint = domaine_VEF.premiere_face_int();
2047 int nb_faces = domaine_VEF.nb_faces_tot();
2048 int nb_arete = domaine.nb_aretes_tot();
2049 int elem1, elem2, face, ok;
2050 IntLists voisins(nb_arete);
2051 DoubleLists coeffs(nb_arete);
2052 ArrOfInt sommets(dimension+2);
2053 ArrOfInt face_opp1(dimension+2);
2054 ArrOfInt face_opp2(dimension+2);
2055 // Boundary faces:
2056 for(int i=0; i<les_cl.size(); i++)
2057 {
2058 const Cond_lim& la_cl = les_cl[i];
2059 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
2060 int nb_faces_bord_tot = le_bord.nb_faces_tot();
2061 for(int ind_face=0; ind_face<nb_faces_bord_tot; ind_face++)
2062 {
2063 ok=okface(ind_face, face, la_cl);
2064 if (ok==-1) break;
2065 elem1=face_voisins(face, 0);
2066 elem2=face_voisins(face, 1);
2067 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets,
2068 face_opp1, face_opp2);
2069 if(ok==3)
2070 {
2071 contribuer_matrice_NeumannPaPa(domaine_VEF, elem1, sommets, voisins, coeffs);
2072 }
2073 else if(ok==4)
2074 {
2075 contribuer_matrice_SymetriePaPa(domaine_VEF, elem1, sommets, voisins, coeffs);
2076 }
2077 else
2078 contribuer_matricePaPa(domaine_VEF, elem1, elem2, sommets, voisins, coeffs);
2079 }
2080 }
2081 face_associee=-1;
2082 for(face=nint; face<nb_faces; face++)
2083 {
2084 elem1=face_voisins(face, 0);
2085 elem2=face_voisins(face, 1);
2086 if (!domaine_VEF.est_une_face_virt_bord(face)) // Process internal faces only
2087 {
2088 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
2089 contribuer_matricePaPa(domaine_VEF, elem1, elem2, sommets, voisins, coeffs);
2090 }
2091 }
2092
2093 DoubleVect diag(nb_arete);
2094 diag = 1;
2095 matrice.typer("Matrice_Bloc");
2096 Matrice_Bloc& matrice_bloc=ref_cast(Matrice_Bloc, matrice.valeur());
2097 matrice_bloc.remplir(voisins, coeffs, diag, domaine.nb_aretes(), domaine.nb_aretes_tot());
2098 Cerr << "Assemblage Pa OK" << finl;
2099}
2100
2101void updatePaPa(const Domaine_dis_base& z,
2102 const Domaine_Cl_dis_base& zcl,
2103 Matrice& matrice,
2104 const DoubleTab& inverse_quantitee_entrelacee)
2105{
2106 int dimension=Objet_U::dimension;
2107 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, z);
2108 const Domaine& domaine=domaine_VEF.domaine();
2109 const Domaine_Cl_VEF& domaine_Cl_VEF=ref_cast(Domaine_Cl_VEF, zcl);
2110 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
2111 const IntTab& face_voisins = domaine_VEF.face_voisins();
2112 int nint = domaine_VEF.premiere_face_int();
2113 int nb_faces = domaine_VEF.nb_faces_tot();
2114 int nb_arete = domaine.nb_aretes();
2115 int elem1, elem2, face, ok;
2116 ArrOfInt sommets(dimension+2);
2117 ArrOfInt face_opp1(dimension+2);
2118 ArrOfInt face_opp2(dimension+2);
2119 Matrice_Bloc& A=ref_cast(Matrice_Bloc, matrice.valeur());
2120 Matrice_Morse_Sym& ARR=ref_cast(Matrice_Morse_Sym, A.get_bloc(0,0).valeur());
2121 Matrice_Morse& ARV=ref_cast(Matrice_Morse, A.get_bloc(0,1).valeur());
2122 Matrice_Morse& AVR=ref_cast(Matrice_Morse, A.get_bloc(1,0).valeur());
2123 Matrice_Morse& AVV=ref_cast(Matrice_Morse, A.get_bloc(1,1).valeur());
2124 // Boundary faces:
2125 for (auto &itr : les_cl)
2126 {
2127 const Cond_lim& la_cl = itr;
2128 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
2129 int nb_faces_bord_tot = le_bord.nb_faces_tot();
2130 for (int ind_face = 0; ind_face < nb_faces_bord_tot; ind_face++)
2131 {
2132 ok = okface(ind_face, face, la_cl);
2133 if (ok == -1)
2134 break;
2135 elem1 = face_voisins(face, 0);
2136 elem2 = face_voisins(face, 1);
2137 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
2138 if (ok == 3)
2139 update_matrice_NeumannPaPa(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, sommets, face_opp1, ARR, ARV, AVR, AVV);
2140 else if (ok == 4)
2141 update_matrice_SymetriePaPa(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, sommets, face_opp1, ARR, ARV, AVR, AVV);
2142 else
2143 update_matricePaPa(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, elem2, sommets, face_opp1, face_opp2, ARR, ARV, AVR, AVV);
2144 }
2145 }
2146 face_associee=-1;
2147 for (face = nint; face < nb_faces; face++)
2148 {
2149 elem1 = face_voisins(face, 0);
2150 elem2 = face_voisins(face, 1);
2151 if (!domaine_VEF.est_une_face_virt_bord(face)) // Process internal faces only
2152 {
2153 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
2154 update_matricePaPa(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, elem2, sommets, face_opp1, face_opp2, ARR, ARV, AVR, AVV);
2155 }
2156 }
2157 for(int i=0; i<nb_arete; i++)
2158 if(ARR(i,i)==0)
2159 {
2160 // Not displayed because too frequent on large meshes
2161 //Cerr << "Modifying orphan edge row " << i << finl;
2162 ARR(i,i)=1;
2163 }
2164
2165 Cerr << "Update Pa OK" << finl;
2166}
2167
2168void assemblerP0Pa(const Domaine_dis_base& z,
2169 const Domaine_Cl_dis_base& zcl,
2170 Matrice& matrice,
2171 const DoubleTab& inverse_quantitee_entrelacee)
2172{
2173 int dimension=Objet_U::dimension;
2174 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, z);
2175 const Domaine& domaine=domaine_VEF.domaine();
2176 const Domaine_Cl_VEF& domaine_Cl_VEF=ref_cast(Domaine_Cl_VEF, zcl);
2177 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
2178 const IntTab& face_voisins = domaine_VEF.face_voisins();
2179 int nint = domaine_VEF.premiere_face_int();
2180 int nb_faces = domaine_VEF.nb_faces_tot();
2181 int nb_elem = domaine.nb_elem_tot();
2182 int elem1, elem2, face, ok;
2183 IntLists voisins(nb_elem);
2184 DoubleLists coeffs(nb_elem);
2185 ArrOfInt sommets(dimension+2);
2186 ArrOfInt face_opp1(dimension+2);
2187 ArrOfInt face_opp2(dimension+2);
2188 // Boundary faces:
2189 for(int i=0; i<les_cl.size(); i++)
2190 {
2191 const Cond_lim& la_cl = les_cl[i];
2192 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
2193 int nb_faces_bord_tot = le_bord.nb_faces_tot();
2194 for(int ind_face=0; ind_face<nb_faces_bord_tot; ind_face++)
2195 {
2196 ok=okface(ind_face, face, la_cl);
2197 if (ok==-1) break;
2198 elem1=face_voisins(face, 0);
2199 elem2=face_voisins(face, 1);
2200 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets,
2201 face_opp1, face_opp2);
2202 if(ok==3)
2203 {
2204 contribuer_matrice_NeumannP0Pa(domaine_VEF, elem1, sommets, voisins, coeffs);
2205 }
2206 else if(ok==4)
2207 {
2208 ; // NOTHING
2209 }
2210 else
2211 contribuer_matriceP0Pa(domaine_VEF, elem1, elem2, sommets, voisins, coeffs);
2212 }
2213 }
2214 face_associee=-1;
2215 for(face=nint; face<nb_faces; face++)
2216 {
2217 elem1=face_voisins(face, 0);
2218 elem2=face_voisins(face, 1);
2219 if (!domaine_VEF.est_une_face_virt_bord(face)) // Process internal faces only
2220 {
2221 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets,
2222 face_opp1, face_opp2);
2223 contribuer_matriceP0Pa(domaine_VEF, elem1, elem2, sommets, voisins, coeffs);
2224 }
2225 }
2226 matrice.typer("Matrice_Bloc");
2227 Matrice_Bloc& matrice_bloc=ref_cast(Matrice_Bloc, matrice.valeur());
2228 matrice_bloc.remplir(voisins, coeffs, domaine.nb_elem(), domaine.nb_elem_tot(), domaine.nb_aretes(), domaine.nb_aretes_tot());
2229 Cerr << "Assemblage P0Pa OK" << finl;
2230}
2231
2232void updateP0Pa(const Domaine_dis_base& z,
2233 const Domaine_Cl_dis_base& zcl,
2234 Matrice& matrice,
2235 const DoubleTab& inverse_quantitee_entrelacee)
2236{
2237 int dimension=Objet_U::dimension;
2238 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, z);
2239 const Domaine_Cl_VEF& domaine_Cl_VEF=ref_cast(Domaine_Cl_VEF, zcl);
2240 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
2241 const IntTab& face_voisins = domaine_VEF.face_voisins();
2242 int nint = domaine_VEF.premiere_face_int();
2243 int nb_faces = domaine_VEF.nb_faces_tot();
2244 int elem1, elem2, face, ok;
2245 ArrOfInt sommets(dimension+2);
2246 ArrOfInt face_opp1(dimension+2);
2247 ArrOfInt face_opp2(dimension+2);
2248 Matrice_Bloc& A=ref_cast(Matrice_Bloc, matrice.valeur());
2249 Matrice_Morse& ARR=ref_cast(Matrice_Morse, A.get_bloc(0,0).valeur());
2250 Matrice_Morse& ARV=ref_cast(Matrice_Morse, A.get_bloc(0,1).valeur());
2251 Matrice_Morse& AVR=ref_cast(Matrice_Morse, A.get_bloc(1,0).valeur());
2252 Matrice_Morse& AVV=ref_cast(Matrice_Morse, A.get_bloc(1,1).valeur());
2253 // Boundary faces:
2254 for (auto &itr : les_cl)
2255 {
2256 const Cond_lim& la_cl = itr;
2257 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
2258 int nb_faces_bord_tot = le_bord.nb_faces_tot();
2259 for (int ind_face = 0; ind_face < nb_faces_bord_tot; ind_face++)
2260 {
2261 ok = okface(ind_face, face, la_cl);
2262 if (ok == -1)
2263 break;
2264 elem1 = face_voisins(face, 0);
2265 elem2 = face_voisins(face, 1);
2266 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
2267 if (ok == 3)
2268 update_matrice_NeumannP0Pa(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, sommets, face_opp1, ARR, ARV, AVR, AVV);
2269 else if (ok == 4) { /* Do nothing */ }
2270 else
2271 update_matriceP0Pa(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, elem2, sommets, face_opp1, face_opp2, ARR, ARV, AVR, AVV);
2272 }
2273 }
2274 face_associee = -1;
2275 for (face = nint; face < nb_faces; face++)
2276 {
2277 elem1 = face_voisins(face, 0);
2278 elem2 = face_voisins(face, 1);
2279 if (!domaine_VEF.est_une_face_virt_bord(face)) // Process internal faces only
2280 {
2281 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
2282 update_matriceP0Pa(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, elem2, sommets, face_opp1, face_opp2, ARR, ARV, AVR, AVV);
2283 }
2284 }
2285 Cerr << "Update P0Pa OK" << finl;
2286}
2287
2288void assemblerP1Pa(const Domaine_dis_base& z,
2289 const Domaine_Cl_dis_base& zcl,
2290 Matrice& matrice,
2291 const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som)
2292{
2293 int dimension=Objet_U::dimension;
2294 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, z);
2295 const Domaine& domaine=domaine_VEF.domaine();
2296 const Domaine_Cl_VEF& domaine_Cl_VEF=ref_cast(Domaine_Cl_VEF, zcl);
2297 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
2298 const IntTab& face_voisins = domaine_VEF.face_voisins();
2299 int nint = domaine_VEF.premiere_face_int();
2300 int nb_faces = domaine_VEF.nb_faces_tot();
2301 int elem1, elem2, face, ok;
2302 int nb_som = domaine.nb_som_tot();
2303 IntLists voisins(nb_som);
2304 DoubleLists coeffs(nb_som);
2305 ArrOfInt sommets(dimension+2);
2306 ArrOfInt face_opp1(dimension+2);
2307 ArrOfInt face_opp2(dimension+2);
2308 // Boundary faces:
2309 for(int i=0; i<les_cl.size(); i++)
2310 {
2311 const Cond_lim& la_cl = les_cl[i];
2312 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
2313 int nb_faces_bord_tot = le_bord.nb_faces_tot();
2314 for(int ind_face=0; ind_face<nb_faces_bord_tot; ind_face++)
2315 {
2316 ok=okface(ind_face, face, la_cl);
2317 if (ok==-1) break;
2318 elem1=face_voisins(face, 0);
2319 elem2=face_voisins(face, 1);
2320 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets,
2321 face_opp1, face_opp2);
2322 if(ok==3)
2323 {
2324 contribuer_matrice_NeumannP1Pa(domaine_VEF, elem1, sommets, voisins, coeffs);
2325 }
2326 else if(ok==4)
2327 {
2328 contribuer_matrice_SymetrieP1Pa(domaine_VEF, elem1, sommets, voisins, coeffs);
2329 }
2330 else
2331 contribuer_matriceP1Pa(domaine_VEF, elem1, elem2, sommets, voisins, coeffs);
2332 }
2333 }
2334 face_associee=-1;
2335 for(face=nint; face<nb_faces; face++)
2336 {
2337 elem1=face_voisins(face, 0);
2338 elem2=face_voisins(face, 1);
2339 if (!domaine_VEF.est_une_face_virt_bord(face)) // Process internal faces only
2340 {
2341 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets,
2342 face_opp1, face_opp2);
2343 contribuer_matriceP1Pa(domaine_VEF, elem1, elem2, sommets, voisins, coeffs);
2344 }
2345 }
2346
2347 matrice.typer("Matrice_Bloc");
2348 Matrice_Bloc& matrice_bloc=ref_cast(Matrice_Bloc, matrice.valeur());
2349 matrice_bloc.remplir(voisins, coeffs, domaine.nb_som(), domaine.nb_som_tot(), domaine.nb_aretes(), domaine.nb_aretes_tot());
2350 Cerr << "Assemblage P1Pa OK" << finl;
2351}
2352
2353void updateP1Pa(const Domaine_dis_base& z,
2354 const Domaine_Cl_dis_base& zcl,
2355 Matrice& matrice,
2356 const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som)
2357{
2358 int dimension=Objet_U::dimension;
2359 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, z);
2360 //const Domaine& domaine=domaine_VEF.domaine();
2361 const Domaine_Cl_VEF& domaine_Cl_VEF=ref_cast(Domaine_Cl_VEF, zcl);
2362 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
2363 const IntTab& face_voisins = domaine_VEF.face_voisins();
2364 int nint = domaine_VEF.premiere_face_int();
2365 int nb_faces = domaine_VEF.nb_faces_tot();
2366 int elem1, elem2, face, ok;
2367 //int nb_som = domaine.nb_som_tot();
2368 //IntLists voisins(nb_som);
2369 //DoubleLists coeffs(nb_som);
2370 ArrOfInt sommets(dimension+2);
2371 ArrOfInt face_opp1(dimension+2);
2372 ArrOfInt face_opp2(dimension+2);
2373 Matrice_Bloc& A=ref_cast(Matrice_Bloc, matrice.valeur());
2374 Matrice_Morse& ARR=ref_cast(Matrice_Morse, A.get_bloc(0,0).valeur());
2375 Matrice_Morse& ARV=ref_cast(Matrice_Morse, A.get_bloc(0,1).valeur());
2376 Matrice_Morse& AVR=ref_cast(Matrice_Morse, A.get_bloc(1,0).valeur());
2377 Matrice_Morse& AVV=ref_cast(Matrice_Morse, A.get_bloc(1,1).valeur());
2378 // Boundary faces:
2379 for (auto &itr : les_cl)
2380 {
2381 const Cond_lim& la_cl = itr;
2382 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
2383 int nb_faces_bord_tot = le_bord.nb_faces_tot();
2384 for (int ind_face = 0; ind_face < nb_faces_bord_tot; ind_face++)
2385 {
2386 ok = okface(ind_face, face, la_cl);
2387 if (ok == -1)
2388 break;
2389 elem1 = face_voisins(face, 0);
2390 elem2 = face_voisins(face, 1);
2391 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
2392 if (ok == 3)
2393 update_matrice_NeumannP1Pa(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, sommets, face_opp1, coef_som, ARR, ARV, AVR, AVV);
2394 else if (ok == 4)
2395 update_matrice_SymetrieP1Pa(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, sommets, face_opp1, coef_som, ARR, ARV, AVR, AVV);
2396 else
2397 update_matriceP1Pa(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, elem2, sommets, face_opp1, face_opp2, coef_som, ARR, ARV, AVR, AVV);
2398 }
2399 }
2400 face_associee=-1;
2401 for (face = nint; face < nb_faces; face++)
2402 {
2403 elem1 = face_voisins(face, 0);
2404 elem2 = face_voisins(face, 1);
2405 if (!domaine_VEF.est_une_face_virt_bord(face)) // Process internal faces only
2406 {
2407 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
2408 update_matriceP1Pa(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, elem2, sommets, face_opp1, face_opp2, coef_som, ARR, ARV, AVR, AVV);
2409 }
2410 }
2411
2412 Cerr << "Update P1Pa OK" << finl;
2413}
2414
2415void assemblerP0P1(const Domaine_dis_base& z,
2416 const Domaine_Cl_dis_base& zcl,
2417 Matrice& matrice,
2418 const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som)
2419{
2420 int dimension=Objet_U::dimension;
2421 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, z);
2422 const Domaine& domaine=domaine_VEF.domaine();
2423 const Domaine_Cl_VEF& domaine_Cl_VEF=ref_cast(Domaine_Cl_VEF, zcl);
2424 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
2425 const IntTab& face_voisins = domaine_VEF.face_voisins();
2426 int nint = domaine_VEF.premiere_face_int();
2427 int nb_faces = domaine_VEF.nb_faces_tot();
2428 int nb_elem = domaine.nb_elem_tot();
2429 int elem1, elem2, face, ok;
2430 IntLists voisins(nb_elem);
2431 DoubleLists coeffs(nb_elem);
2432 ArrOfInt sommets(dimension+2);
2433 ArrOfInt face_opp1(dimension+2);
2434 ArrOfInt face_opp2(dimension+2);
2435 // Boundary faces:
2436 for(int i=0; i<les_cl.size(); i++)
2437 {
2438 const Cond_lim& la_cl = les_cl[i];
2439 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
2440 int nb_faces_bord_tot = le_bord.nb_faces_tot();
2441 for(int ind_face=0; ind_face<nb_faces_bord_tot; ind_face++)
2442 {
2443 ok=okface(ind_face, face, la_cl);
2444 if (ok==-1) break;
2445 elem1=face_voisins(face, 0);
2446 elem2=face_voisins(face, 1);
2447 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
2448 sort(sommets, face_opp1, face_opp2);
2449 if(ok==3)
2450 {
2451 contribuer_matrice_NeumannP0P1(elem1, sommets, voisins, coeffs);
2452 }
2453 else if(ok==4)
2454 {
2455 ;// NOTHING
2456 }
2457 else
2458 contribuer_matriceP0P1(elem1, elem2, sommets, voisins, coeffs);
2459 }
2460 }
2461 face_associee=-1;
2462 for(face=nint; face<nb_faces; face++)
2463 {
2464 elem1=face_voisins(face, 0);
2465 elem2=face_voisins(face, 1);
2466 if (!domaine_VEF.est_une_face_virt_bord(face)) // Process internal faces only
2467 {
2468 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
2469 sort(sommets, face_opp1, face_opp2);
2470 contribuer_matriceP0P1(elem1, elem2, sommets, voisins, coeffs);
2471 }
2472 }
2473
2474 matrice.typer("Matrice_Bloc");
2475 Matrice_Bloc& matrice_bloc=ref_cast(Matrice_Bloc, matrice.valeur());
2476 matrice_bloc.remplir(voisins, coeffs, domaine.nb_elem(), domaine.nb_elem_tot(), domaine.nb_som(), domaine.nb_som_tot());
2477 Cerr << "Assemblage POP1 OK" << finl;
2478}
2479
2480void updateP0P1(const Domaine_dis_base& z,
2481 const Domaine_Cl_dis_base& zcl,
2482 Matrice& matrice,
2483 const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som)
2484{
2485 int dimension=Objet_U::dimension;
2486 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, z);
2487 const Domaine_Cl_VEF& domaine_Cl_VEF=ref_cast(Domaine_Cl_VEF, zcl);
2488 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
2489 const IntTab& face_voisins = domaine_VEF.face_voisins();
2490 int nint = domaine_VEF.premiere_face_int();
2491 int nb_faces = domaine_VEF.nb_faces_tot();
2492 int elem1, elem2, face, ok;
2493 ArrOfInt sommets(dimension+2);
2494 ArrOfInt face_opp1(dimension+2);
2495 ArrOfInt face_opp2(dimension+2);
2496 Matrice_Bloc& A=ref_cast(Matrice_Bloc, matrice.valeur());
2497 Matrice_Morse& ARR=ref_cast(Matrice_Morse, A.get_bloc(0,0).valeur());
2498 Matrice_Morse& ARV=ref_cast(Matrice_Morse, A.get_bloc(0,1).valeur());
2499 Matrice_Morse& AVR=ref_cast(Matrice_Morse, A.get_bloc(1,0).valeur());
2500 Matrice_Morse& AVV=ref_cast(Matrice_Morse, A.get_bloc(1,1).valeur());
2501 // Boundary faces:
2502 for (auto &itr : les_cl)
2503 {
2504 const Cond_lim& la_cl = itr;
2505 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
2506 int nb_faces_bord_tot = le_bord.nb_faces_tot();
2507 for (int ind_face = 0; ind_face < nb_faces_bord_tot; ind_face++)
2508 {
2509 ok = okface(ind_face, face, la_cl);
2510 if (ok == -1)
2511 break;
2512 elem1 = face_voisins(face, 0);
2513 elem2 = face_voisins(face, 1);
2514 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets, face_opp1, face_opp2);
2515 sort(sommets, face_opp1, face_opp2);
2516 if (ok == 3)
2517 update_matrice_NeumannP0P1(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, sommets, face_opp1, coef_som, ARR, ARV, AVR, AVV);
2518 else if (ok == 4) { /* Do nothing */ }
2519 else
2520 update_matriceP0P1(domaine_VEF, inverse_quantitee_entrelacee, face, elem1, elem2, sommets, face_opp1, face_opp2, coef_som, ARR, ARV, AVR, AVV);
2521 }
2522 }
2523 face_associee=-1;
2524 for(face=nint; face<nb_faces; face++)
2525 {
2526 elem1=face_voisins(face, 0);
2527 elem2=face_voisins(face, 1);
2528 if (!domaine_VEF.est_une_face_virt_bord(face)) // Process internal faces only
2529 {
2530 remplir_sommets(domaine_VEF, face, elem1, elem2, sommets,
2531 face_opp1, face_opp2);
2532 sort(sommets, face_opp1, face_opp2);
2533 update_matriceP0P1(domaine_VEF,
2534 inverse_quantitee_entrelacee,
2535 face, elem1, elem2, sommets,
2536 face_opp1, face_opp2, coef_som,
2537 ARR,ARV,AVR,AVV);
2538 }
2539 }
2540 Cerr << "Update POP1 OK" << finl;
2541}
2542
2543
Assembleur_P_VEFPreP1B
Definition Assembleur_P_VEFPreP1B.h:26
Assembleur_P_VEFPreP1B::equation
const Equation_base & equation() const
Definition Assembleur_P_VEFPreP1B.h:41
Assembleur_P_VEF
Definition Assembleur_P_VEF.h:28
Assembleur_P_VEF::remplir
int remplir(Matrice &, const DoubleTab &)
Definition Assembleur_P_VEF.cpp:138
Assembleur_P_VEF::associer_domaine_cl_dis_base
void associer_domaine_cl_dis_base(const Domaine_Cl_dis_base &) override
Definition Assembleur_P_VEF.cpp:698
Assembleur_P_VEF::associer_domaine_dis_base
void associer_domaine_dis_base(const Domaine_dis_base &) override
Definition Assembleur_P_VEF.cpp:693
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Returns the array of field values at the current time.
Definition Champ_Inc_base.h:77
Cond_lim_base
class Cond_lim_base Base class for the hierarchy of classes that represent the different boundary con...
Definition Cond_lim_base.h:40
Cond_lim_base::frontiere_dis
virtual Frontiere_dis_base & frontiere_dis()
Returns the discretized boundary to which the boundary conditions apply.
Definition Cond_lim_base.h:101
Cond_lim
class Cond_lim Generic class used to represent any class
Definition Cond_lim.h:31
Conds_lim
class Conds_lim This class represents a vector of boundary conditions.
Definition Conds_lim.h:32
Dirichlet_homogene
Classe Dirichlet_homogene This class is the base class of the hierarchy of homogeneous Dirichlet-type...
Definition Dirichlet_homogene.h:31
Dirichlet
Dirichlet This class is the base class of the hierarchy of Dirichlet-type boundary conditions.
Definition Dirichlet.h:31
Domaine_32_64::aretes_som
const IntTab_t & aretes_som() const
returns the connectivity array edges/vertices.
Definition Domaine.h:156
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_elem
int_t nb_elem() const
Definition Domaine.h:131
Domaine_32_64::nb_aretes
int_t nb_aretes() const
Returns the number of real edges.
Definition Domaine.h:143
Domaine_32_64::elem_aretes
int_t elem_aretes(int_t i, int j) const
returns the number of the j-th edge of the i-th element.
Definition Domaine.h:154
Domaine_32_64::nb_som_tot
int_t nb_som_tot() const
Returns the total number of vertices of the domain i.e. the number of real and virtual vertices on th...
Definition Domaine.h:123
Domaine_32_64::nb_som
int_t nb_som() const
Returns the number of vertices of the domain.
Definition Domaine.h:121
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:38
Domaine_Cl_dis_base
class Domaine_Cl_dis_base Domaine_Cl_dis_base objects represent discretized boundary conditions
Definition Domaine_Cl_dis_base.h:37
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Returns the i-th boundary condition.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:53
Domaine_VEF::get_ok_arete
const IntVect & get_ok_arete() const
Definition Domaine_VEF.h:98
Domaine_VEF::get_alphaA
int get_alphaA() const
Definition Domaine_VEF.h:93
Domaine_VEF::get_renum_arete_perio
const ArrOfInt & get_renum_arete_perio() const
Definition Domaine_VEF.h:97
Domaine_VEF::get_alphaS
int get_alphaS() const
Definition Domaine_VEF.h:92
Domaine_VEF::get_alphaE
int get_alphaE() const
Definition Domaine_VEF.h:91
Domaine_VF::nb_faces_tot
int nb_faces_tot() const
Returns the total number of 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::face_sommets
int face_sommets(int i, int j) const
Returns the index of the i-th vertex of face num_face.
Definition Domaine_VF.h:582
Domaine_VF::est_une_face_virt_bord
int est_une_face_virt_bord(int) const
Returns 1 if face is a virtual boundary face, 0 otherwise.
Definition Domaine_VF.h:680
Domaine_VF::elem_faces
int elem_faces(int i, int j) const
Returns the index of the i-th face of element num_elem; the face numbering convention is.
Definition Domaine_VF.h:542
Domaine_VF::premiere_face_int
int premiere_face_int() const
A face is internal if and only if it separates two elements.
Definition Domaine_VF.h:463
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
Returns the neighbouring element of num_face in direction i.
Definition Domaine_VF.h:418
Domaine_dis_base
class Domaine_dis_base This class is the base of the hierarchy of discretized domains.
Definition Domaine_dis_base.h:39
Domaine_dis_base::nb_elem
int nb_elem() const
Definition Domaine_dis_base.h:49
Domaine_dis_base::nb_som
int nb_som() const
Definition Domaine_dis_base.h:51
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Equation_base::solv_masse
Solveur_Masse_base & solv_masse()
Returns the mass solver associated with the equation.
Definition Equation_base.h:363
Front_VF
class Front_VF
Definition Front_VF.h:36
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
Matrice_Base::multvect
virtual DoubleVect & multvect(const DoubleVect &, DoubleVect &) const
Multiplication of a vector by the matrix.
Definition Matrice_Base.h:103
Matrice_Bloc_Sym
Definition Matrice_Bloc_Sym.h:62
Matrice_Bloc_Sym::imprimer_formatte
Sortie & imprimer_formatte(Sortie &s) const override
Definition Matrice_Bloc_Sym.cpp:57
Matrice_Bloc
Definition Matrice_Bloc.h:51
Matrice_Bloc::get_bloc
virtual const Matrice & get_bloc(int i, int j) const
Definition Matrice_Bloc.cpp:601
Matrice_Bloc::remplir
void remplir(const IntLists &voisins, const DoubleLists &valeurs, const DoubleVect &terme_diag, const int i, const int n)
Definition Matrice_Bloc.cpp:693
Matrice_Morse_Sym
Matrice_Morse_Sym class - Represents a sparse symmetric matrix M stored in Morse format.
Definition Matrice_Morse_Sym.h:34
Matrice_Morse
Matrice_Morse class - Represents a (sparse) matrix M, not necessarily square,.
Definition Matrice_Morse.h:50
Matrice
Matrice class - Generic class in the matrix hierarchy.
Definition Matrice.h:34
Navier_Stokes_std
Navier_Stokes_std This class carries the terms of the momentum equation.
Definition Navier_Stokes_std.h:55
Navier_Stokes_std::inconnue
const Champ_Inc_base & inconnue() const override
Returns the velocity (unknown field of the equation) (const version).
Definition Navier_Stokes_std.cpp:597
Navier_Stokes_std::operateur_divergence
Operateur_Div & operateur_divergence()
Returns the divergence operator associated with the equation.
Definition Navier_Stokes_std.cpp:546
Navier_Stokes_std::operateur_gradient
Operateur_Grad & operateur_gradient()
Returns the gradient operator associated with the equation.
Definition Navier_Stokes_std.cpp:566
Navier_Stokes_std::pression
Champ_Inc_base & pression()
Definition Navier_Stokes_std.h:122
Neumann_sortie_libre
Neumann_sortie_libre This class represents an open boundary without imposed velocity.
Definition Neumann_sortie_libre.h:34
Objet_U::dimension
static int dimension
Definition Objet_U.h:94
Op_Div_VEFP1B_Elem
class Op_Div_VEFP1B_Elem
Definition Op_Div_VEFP1B_Elem.h:36
Op_Grad_VEF_P1B_Face
class Op_Grad_VEF_P1B_Face
Definition Op_Grad_VEF_P1B_Face.h:33
Operateur_Div
Operateur_Div Generic class of the hierarchy of operators computing the divergence.
Definition Operateur_Div.h:31
Operateur_Grad
Classe Operateur_Grad Generic class of the hierarchy of operators computing the gradient.
Definition Operateur_Grad.h:31
Operateur_base::calculer
virtual DoubleTab & calculer(const DoubleTab &, DoubleTab &) const
Definition Operateur_base.cpp:242
Periodique
class Periodique This class represents a periodic boundary condition.
Definition Periodique.h:31
Periodique::face_associee
int face_associee(int i) const
Definition Periodique.h:35
Process::mp_sum_for_each
static void mp_sum_for_each(T &arg1, T &arg2)
C++14 compatible mp_sum_for_each: combine multiple mp_sum calls into one collective operation Usage: ...
Definition Process.cpp:208
Process::nproc
static int nproc()
Returns the number of processors in the current group. See Comm_Group::nproc() and PE_Groups::current...
Definition Process.cpp:102
Process::me
static int me()
Returns the rank of the local processor in the current communication group. See Comm_Group::rank() an...
Definition Process.cpp:122
Process::exit
static void exit(int exit_code=-1)
Exit routine for TRUST within a Kokkos region.
Definition Process.cpp:466
Solveur_Masse_base
Solveur_Masse_base Represents the mass matrix of an equation.
Definition Solveur_Masse_base.h:40
Solveur_Masse_base::appliquer
virtual DoubleTab & appliquer(DoubleTab &) const
Returns appliquer_impl(x/temporal_coefficient) if a temporal coefficient is set, otherwise returns ap...
Definition Solveur_Masse_base.cpp:91
Symetrie
Symetrie On symmetry faces, the following properties hold:
Definition Symetrie.h:37
TRUSTArray::size_array
_SIZE_ size_array() const
Definition TRUSTArray.tpp:187
TRUSTTab::ecrit
void ecrit(Sortie &) const override
Definition TRUSTTab.tpp:688
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133
TRUSTVect::size
_SIZE_ size() const
Definition TRUSTVect.tpp:45