LORENE
poisson_vect_frontiere.C
1/*
2 * Copyright (c) 2000-2001 Philippe Grandclement
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23char poisson_vect_frontiere_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_vect_frontiere.C,v 1.9 2014/10/13 08:53:30 j_novak Exp $" ;
24
25/*
26 * $Id: poisson_vect_frontiere.C,v 1.9 2014/10/13 08:53:30 j_novak Exp $
27 * $Log: poisson_vect_frontiere.C,v $
28 * Revision 1.9 2014/10/13 08:53:30 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.8 2014/10/06 15:16:09 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.7 2005/03/11 11:20:26 f_limousin
35 * Minor modif
36 *
37 * Revision 1.6 2005/02/22 18:00:32 f_limousin
38 * Correction of an error in the function poisson_vect_binaire(...).
39 * Confusion between cartesian and spherical triad for the solution.
40 *
41 * Revision 1.5 2005/02/08 10:07:07 f_limousin
42 * Implementation of poisson_vect_binaire(...) with Vectors (instead of
43 * Tenseur) in argument.
44 *
45 * Revision 1.4 2004/09/28 16:00:15 f_limousin
46 * Add function poisson_vect_boundary which is the same as
47 * poisson_vect_frontiere but for the new classes Tensor and Scalar.
48 *
49 * Revision 1.3 2003/10/03 15:58:50 j_novak
50 * Cleaning of some headers
51 *
52 * Revision 1.2 2003/02/13 16:40:25 p_grandclement
53 * Addition of various things for the Bin_ns_bh project, non of them being
54 * completely tested
55 *
56 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
57 * LORENE
58 *
59 * Revision 2.2 2000/10/26 09:08:06 phil
60 * *** empty log message ***
61 *
62 * Revision 2.1 2000/10/26 09:01:18 phil
63 * *** empty log message ***
64 *
65 * Revision 2.0 2000/10/19 09:36:36 phil
66 * *** empty log message ***
67 *
68 *
69 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_vect_frontiere.C,v 1.9 2014/10/13 08:53:30 j_novak Exp $
70 *
71 */
72
73// Header C :
74#include <cstdlib>
75#include <cmath>
76
77// Headers Lorene :
78#include "proto.h"
79#include "tenseur.h"
80#include "tensor.h"
81#include "metric.h"
82
83 // USING OOhara
84namespace Lorene {
85void poisson_vect_frontiere (double lambda, const Tenseur& source, Tenseur& shift,
86 const Valeur& lim_x, const Valeur& lim_y, const Valeur& lim_z,
87 int num_front, double precision, int itermax) {
88
89 // METTRE TOUT PLEIN D'ASSERT
90
91 // Confort
92 int nt = lim_x.get_mg()->get_nt(num_front+1) ;
93 int np = lim_x.get_mg()->get_np(num_front+1) ;
94 int nz = lim_x.get_mg()->get_nzone() ;
95
96 if (shift.get_etat() == ETATZERO) {
97 shift.set_etat_qcq() ;
98 for (int i=0 ; i<3 ; i++)
99 shift.set(i).annule_hard() ;
100 shift.set_std_base() ;
101 }
102
103 Tenseur so (source) ;
104
105 // La source scalaire :
106 Tenseur cop_so (so) ;
107 cop_so.dec2_dzpuis() ;
108 cop_so.dec2_dzpuis() ;
109
110 Tenseur scal (*so.get_mp()) ;
111 scal.set_etat_qcq() ;
112
113 Cmp source_scal (contract(cop_so.gradient(), 0, 1)()/(lambda+1)) ;
114 source_scal.inc2_dzpuis() ;
115 if (source_scal.get_etat()== ETATZERO) {
116 source_scal.annule_hard() ;
117 source_scal.std_base_scal() ;
118 source_scal.set_dzpuis(4) ;
119 }
120
121 Tenseur copie_so (so) ;
122 copie_so.dec_dzpuis() ;
123
124 Tenseur source_vect (*so.get_mp(), 1, CON, *source.get_triad()) ;
125 Tenseur auxi (*so.get_mp(), 1, COV, *source.get_triad()) ;
126 Cmp grad_shift (source_scal.get_mp()) ;
127
128 // La condition sur la derivee du scalaire :
129 Valeur lim_scal (lim_x.get_mg()) ;
130 Tenseur shift_old (*shift.get_mp(), 1, CON, shift.get_mp()->get_bvect_cart()) ;
131
132 int conte = 0 ;
133 int indic = 1 ;
134
135 while (indic ==1) {
136
137 shift_old = shift ;
138
139 grad_shift = contract(shift.gradient(), 0, 1)() ;
140 grad_shift.dec2_dzpuis() ;
141 grad_shift.va.coef_i() ;
142
143
144
145 lim_scal = 1 ; // Permet d'affecter les trucs qui vont bien !
146 for (int k=0 ; k<np ; k++)
147 for (int j=0 ; j<nt ; j++)
148 lim_scal.set(num_front, k, j, 0) =
149 grad_shift.va (num_front+1, k, j, 0) ;
150 lim_scal.std_base_scal() ;
151
152 // On resout la scalaire :
153 scal.set() = source_scal.poisson_dirichlet (lim_scal, num_front) ;
154
155 // La source vectorielle :
156 source_vect.set_etat_qcq() ;
157 auxi = scal.gradient() ;
158 auxi.inc_dzpuis() ;
159 for (int i=0 ; i<3 ; i++)
160 source_vect.set(i) = copie_so(i) - lambda * auxi(i) ;
161
162 indic = 0;
163 for (int i=0 ; i<3 ; i++)
164 if (source_vect(i).get_etat()==ETATQCQ)
165 indic = 1 ;
166 if (indic==0) {
167 for (int i=0 ; i<3 ; i++)
168 source_vect.set(i).annule_hard() ;
169 source_vect.set_std_base() ;
170 }
171
172 // On resout les equations de poisson sur le shift :
173 shift.set(0) = source_vect(0).poisson_dirichlet (lim_x, num_front) ;
174 shift.set(1) = source_vect(1).poisson_dirichlet (lim_y, num_front) ;
175 shift.set(2) = source_vect(2).poisson_dirichlet (lim_z, num_front) ;
176
177 double erreur = 0 ;
178 for (int i=0 ; i<3 ; i++)
179 if (max(norme(shift(i))) > precision) {
180 Tbl diff (diffrelmax (shift(i), shift_old(i))) ;
181 for (int j=num_front+1 ; j<nz ; j++)
182 if (diff(j)> erreur)
183 erreur = diff(j) ;
184 }
185
186 cout << "Pas " << conte << " : Difference " << erreur << endl ;
187 conte ++ ;
188
189 if ((erreur <precision) || (conte > itermax))
190 indic = -1 ;
191 }
192}
193
194
195
196 // USING OOhara
197void poisson_vect_boundary (double lambda, const Vector& source,Vector& shift,
198 const Valeur& lim_x, const Valeur& lim_y, const Valeur& lim_z,
199 int num_front, double precision, int itermax) {
200
201 // On travaille en composantes cartesiennes
202 assert(source.get_mp().get_bvect_spher() == *(source.get_triad())) ;
203 assert(source.get_mp().get_bvect_spher() == *(shift.get_triad())) ;
204
205
206 // Confort
207 int nt = lim_x.get_mg()->get_nt(num_front+1) ;
208 int np = lim_x.get_mg()->get_np(num_front+1) ;
209 int nz = lim_x.get_mg()->get_nzone() ;
210
211 Metric_flat ff(source.get_mp(), source.get_mp().get_bvect_spher()) ;
212
213 Vector so (source) ;
214
215 // La source scalaire :
216 Vector cop_so (so) ;
217 cop_so.dec_dzpuis(2) ;
218 cop_so.dec_dzpuis(2) ;
219
220 Scalar scal (so.get_mp()) ;
221
222 Scalar source_scal (contract(cop_so.derive_cov(ff), 0, 1)/(lambda+1)) ;
223 source_scal.inc_dzpuis(2) ;
224 if (source_scal.get_etat()== ETATZERO) {
225 source_scal.annule_hard() ;
226 source_scal.std_spectral_base() ;
227 source_scal.set_dzpuis(4) ;
228 }
229
230 Vector copie_so (so) ;
231 copie_so.dec_dzpuis() ;
232
233 Vector source_vect (so.get_mp(), CON, *source.get_triad()) ;
234 Vector auxi (so.get_mp(), COV, *source.get_triad()) ;
235 Scalar grad_shift (source_scal.get_mp()) ;
236
237 // La condition sur la derivee du scalaire :
238 Valeur lim_scal (lim_x.get_mg()) ;
239 Vector shift_old (shift.get_mp(), CON, shift.get_mp().get_bvect_cart()) ;
240
241 int conte = 0 ;
242 int indic = 1 ;
243
244 while (indic ==1) {
245
246 shift_old = shift ;
247
248 grad_shift = contract(shift.derive_cov(ff), 0, 1) ;
249 grad_shift.dec_dzpuis(2) ;
250 grad_shift.set_spectral_va().coef_i() ;
251
252
253 lim_scal = 1 ; // Permet d'affecter les trucs qui vont bien !
254 for (int k=0 ; k<np ; k++)
255 for (int j=0 ; j<nt ; j++)
256 lim_scal.set(num_front, k, j, 0) =
257 grad_shift.get_spectral_va() (num_front+1, k, j, 0) ;
258 lim_scal.std_base_scal() ;
259
260 // On resout la scalaire :
261
262 source_scal.filtre(4) ;
263 scal = source_scal.poisson_dirichlet (lim_scal, num_front) ;
264
265 // La source vectorielle :
266 source_vect.set_etat_qcq() ;
267 auxi = scal.derive_cov(ff) ;
268 auxi.inc_dzpuis() ;
269 for (int i=1 ; i<=3 ; i++)
270 source_vect.set(i) = copie_so(i) - lambda * auxi(i) ;
271
272 indic = 0;
273 for (int i=1 ; i<=3 ; i++)
274 if (source_vect(i).get_etat()==ETATQCQ)
275 indic = 1 ;
276 if (indic==0) {
277 for (int i=1 ; i<=3 ; i++)
278 source_vect.set(i).annule_hard() ;
279 source_vect.std_spectral_base() ;
280 }
281
282 shift.change_triad(source.get_mp().get_bvect_cart()) ;
283 source_vect.change_triad(source.get_mp().get_bvect_cart()) ;
284
285 for (int i=1 ; i<=3 ; i++)
286 source_vect.set(i).filtre(4) ;
287
288 // On resout les equations de poisson sur le shift :
289 shift.set(1) = source_vect(1).poisson_dirichlet (lim_x, num_front) ;
290 shift.set(2) = source_vect(2).poisson_dirichlet (lim_y, num_front) ;
291 shift.set(3) = source_vect(3).poisson_dirichlet (lim_z, num_front) ;
292
293 shift.change_triad(source.get_mp().get_bvect_spher()) ;
294 source_vect.change_triad(source.get_mp().get_bvect_spher()) ;
295
296 double erreur = 0 ;
297 for (int i=1 ; i<=3 ; i++)
298 if (max(norme(shift(i))) > precision) {
299 Tbl diff (diffrelmax (shift(i), shift_old(i))) ;
300 for (int j=num_front+1 ; j<nz ; j++)
301 if (diff(j)> erreur)
302 erreur = diff(j) ;
303 }
304
305 cout << "Pas " << conte << " : Difference " << erreur << endl ;
306 conte ++ ;
307
308 if ((erreur <precision) || (conte > itermax))
309 indic = -1 ;
310 }
311}
312
313
314
315
316void poisson_vect_binaire ( double lambda,
317 const Tenseur& source_un, const Tenseur& source_deux,
318 const Valeur& bound_x_un, const Valeur& bound_y_un,
319 const Valeur& bound_z_un, const Valeur& bound_x_deux,
320 const Valeur& bound_y_deux, const Valeur& bound_z_deux,
321 Tenseur& sol_un, Tenseur& sol_deux, int num_front, double precision) {
322
323 // METTRE DES ASSERT
324 assert (sol_un.get_etat() != ETATNONDEF) ;
325 assert (sol_deux.get_etat() != ETATNONDEF) ;
326
327 // Les bases des deux vecteurs doivent etre alignees ou non alignees :
328
329 assert (sol_un.get_mp() == source_un.get_mp()) ;
330 assert (sol_deux.get_mp() == source_deux.get_mp()) ;
331
332 double orientation_un = sol_un.get_mp()->get_rot_phi() ;
333 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
334
335 double orientation_deux = sol_deux.get_mp()->get_rot_phi() ;
336 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
337
338 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
339
340
341 if (sol_un.get_etat() == ETATZERO) {
342 sol_un.set_etat_qcq() ;
343 for (int i=0 ; i<3 ; i++)
344 sol_un.set(i).annule_hard() ;
345 sol_un.set_std_base() ;
346 }
347
348 if (sol_deux.get_etat() == ETATZERO) {
349 sol_deux.set_etat_qcq() ;
350 for (int i=0 ; i<3 ; i++)
351 sol_deux.set(i).annule_hard() ;
352 sol_deux.set_std_base() ;
353 }
354
355 Valeur limite_x_un (bound_x_un.get_mg()) ;
356 limite_x_un = bound_x_un ;
357 Valeur limite_y_un (bound_y_un.get_mg()) ;
358 limite_y_un = bound_y_un ;
359 Valeur limite_z_un (bound_z_un.get_mg()) ;
360 limite_z_un = bound_z_un ;
361
362 Valeur limite_x_deux (bound_x_deux.get_mg()) ;
363 limite_x_deux = bound_x_deux ;
364 Valeur limite_y_deux (bound_y_deux.get_mg()) ;
365 limite_y_deux = bound_y_deux ;
366 Valeur limite_z_deux (bound_z_deux.get_mg()) ;
367 limite_z_deux = bound_z_deux ;
368
369 Valeur limite_chi_un (bound_x_un.get_mg()) ;
370 limite_chi_un = 0 ;
371 limite_chi_un.std_base_scal() ;
372
373 Valeur limite_chi_deux (bound_x_deux.get_mg()) ;
374 limite_chi_deux = 0 ;
375 limite_chi_deux.std_base_scal() ;
376
377 Mtbl xa_mtbl_un (source_un.get_mp()->get_mg()) ;
378 xa_mtbl_un.set_etat_qcq() ;
379 Mtbl ya_mtbl_un (source_un.get_mp()->get_mg()) ;
380 ya_mtbl_un.set_etat_qcq() ;
381 Mtbl za_mtbl_un (source_un.get_mp()->get_mg()) ;
382 za_mtbl_un.set_etat_qcq() ;
383 Mtbl xa_mtbl_deux (source_deux.get_mp()->get_mg()) ;
384 xa_mtbl_deux.set_etat_qcq() ;
385 Mtbl ya_mtbl_deux (source_deux.get_mp()->get_mg()) ;
386 ya_mtbl_deux.set_etat_qcq() ;
387 Mtbl za_mtbl_deux (source_deux.get_mp()->get_mg()) ;
388 za_mtbl_deux.set_etat_qcq() ;
389
390 xa_mtbl_un = source_un.get_mp()->xa ;
391 ya_mtbl_un = source_un.get_mp()->ya ;
392 za_mtbl_un = source_un.get_mp()->za ;
393
394 xa_mtbl_deux = source_deux.get_mp()->xa ;
395 ya_mtbl_deux = source_deux.get_mp()->ya ;
396 za_mtbl_deux = source_deux.get_mp()->za ;
397
398 double xabs, yabs, zabs ;
399 double air, theta, phi ;
400 double valeur ;
401
402 int nbrep_un = bound_x_un.get_mg()->get_np(num_front) ;
403 int nbret_un = bound_x_un.get_mg()->get_nt(num_front) ;
404 int nbrep_deux = bound_x_deux.get_mg()->get_np(num_front) ;
405 int nbret_deux = bound_x_deux.get_mg()->get_nt(num_front) ;
406 int nz_un = bound_x_un.get_mg()->get_nzone() ;
407 int nz_deux = bound_x_deux.get_mg()->get_nzone() ;
408
409 // La source de l'equation scalaire sur 1
410 Tenseur cop_so_un (source_un) ;
411 cop_so_un.dec2_dzpuis() ;
412 cop_so_un.dec2_dzpuis() ;
413
414 Cmp source_scal_un (contract (cop_so_un.gradient(), 0, 1)()/(lambda+1)) ;
415 if (source_scal_un.get_etat() == ETATZERO) {
416 source_scal_un.annule_hard() ;
417 source_scal_un.std_base_scal() ;
418 }
419 source_scal_un.inc2_dzpuis() ;
420
421 // La source de l'equation scalaire sur 2
422 Tenseur cop_so_deux (source_deux) ;
423 cop_so_deux.dec2_dzpuis() ;
424 cop_so_deux.dec2_dzpuis() ;
425
426 Cmp source_scal_deux (contract (cop_so_deux.gradient(), 0, 1)()/(lambda+1)) ;
427 if (source_scal_deux.get_etat() == ETATZERO) {
428 source_scal_deux.annule_hard() ;
429 source_scal_deux.std_base_scal() ;
430 }
431 source_scal_deux.inc2_dzpuis() ;
432
433 // Les copies :
434 Tenseur copie_so_un (source_un) ;
435 copie_so_un.dec_dzpuis() ;
436
437 Tenseur copie_so_deux (source_deux) ;
438 copie_so_deux.dec_dzpuis() ;
439
440 // ON COMMENCE LA BOUCLE :
441 Tenseur sol_un_old (sol_un) ;
442 Tenseur sol_deux_old (sol_deux) ;
443
444 int indic = 1 ;
445 int conte = 0 ;
446
447 while (indic == 1) {
448
449 // On resout les deux equations scalaires :
450 Tenseur chi_un (source_scal_un.poisson_dirichlet (limite_chi_un, num_front)) ;
451 Tenseur chi_deux (source_scal_deux.poisson_dirichlet (limite_chi_deux, num_front)) ;
452
453 // On calcul les source pour les equation vectorielles :
454 Tenseur source_vect_un (copie_so_un) ;
455 if (source_vect_un.get_etat() == ETATZERO) {
456 source_vect_un.set_etat_qcq() ;
457 for (int i=0 ; i<3 ; i++) {
458 source_vect_un.set(i).annule_hard() ;
459 source_vect_un.set(i).set_dzpuis(3) ;
460 }
461 source_vect_un.set_std_base() ;
462 }
463 Tenseur grad_chi_un (chi_un.gradient()) ;
464 grad_chi_un.inc_dzpuis() ;
465 for (int i=0 ; i<3 ; i++)
466 source_vect_un.set(i) = source_vect_un(i)-lambda*grad_chi_un(i) ;
467
468 Tenseur source_vect_deux (copie_so_deux) ;
469 if (source_vect_deux.get_etat() == ETATZERO) {
470 source_vect_deux.set_etat_qcq() ;
471 for (int i=0 ; i<3 ; i++) {
472 source_vect_deux.set(i).annule_hard() ;
473 source_vect_deux.set(i).set_dzpuis(3) ;
474 }
475 source_vect_deux.set_std_base() ;
476 }
477 Tenseur grad_chi_deux (chi_deux.gradient()) ;
478 grad_chi_deux.inc_dzpuis() ;
479 for (int i=0 ; i<3 ; i++)
480 source_vect_deux.set(i) = source_vect_deux(i)-lambda*grad_chi_deux(i) ;
481
482
483 sol_un_old = sol_un ;
484 sol_deux_old = sol_deux ;
485
486 // On resout les equation vectorielles :
487 sol_un.set(0) = source_vect_un(0).poisson_dirichlet (limite_x_un, num_front) ;
488 sol_un.set(1) = source_vect_un(1).poisson_dirichlet (limite_y_un, num_front) ;
489 sol_un.set(2) = source_vect_un(2).poisson_dirichlet (limite_z_un, num_front) ;
490 sol_deux.set(0) = source_vect_deux(0).poisson_dirichlet (limite_x_deux, num_front) ;
491 sol_deux.set(1) = source_vect_deux(1).poisson_dirichlet (limite_y_deux, num_front) ;
492 sol_deux.set(2) = source_vect_deux(2).poisson_dirichlet (limite_z_deux, num_front) ;
493
494
495 // On modifie les Cl sur chi :
496 Cmp div_shift_un (contract(sol_un.gradient(), 0, 1)()) ;
497 div_shift_un.dec2_dzpuis() ;
498 div_shift_un.va.coef_i() ;
499
500 limite_chi_un = 1 ; // Affectation
501 for (int k=0 ; k<nbrep_un ; k++)
502 for (int j=0 ; j<nbret_un ; j++)
503 limite_chi_un.set(num_front, k, j, 0) =
504 div_shift_un.va (num_front+1, k, j, 0) ;
505 limite_chi_un.std_base_scal() ;
506
507 Cmp div_shift_deux (contract(sol_deux.gradient(), 0, 1)()) ;
508 div_shift_deux.dec2_dzpuis() ;
509 div_shift_deux.va.coef_i() ;
510
511 limite_chi_deux = 1 ; // Affectation
512 for (int k=0 ; k<nbrep_deux ; k++)
513 for (int j=0 ; j<nbret_deux ; j++)
514 limite_chi_deux.set(num_front, k, j, 0) =
515 div_shift_deux.va (num_front+1, k, j, 0) ;
516 limite_chi_deux.std_base_scal() ;
517
518
519 // On modifie les Cl sur sol_un :
520 for (int k=0 ; k<nbrep_un ; k++)
521 for (int j=0 ; j<nbret_un ; j++) {
522 xabs = xa_mtbl_un (num_front+1, k, j, 0) ;
523 yabs = ya_mtbl_un (num_front+1, k, j, 0) ;
524 zabs = za_mtbl_un (num_front+1, k, j, 0) ;
525
526 source_deux.get_mp()->convert_absolute
527 (xabs, yabs, zabs, air, theta, phi) ;
528
529 valeur = sol_deux(0).val_point(air, theta, phi) ;
530 limite_x_un.set(num_front, k, j, 0) =
531 bound_x_un(num_front, k, j, 0) - same_orient*valeur ;
532
533 valeur = sol_deux(1).val_point(air, theta, phi) ;
534 limite_y_un.set(num_front, k, j, 0) =
535 bound_y_un(num_front, k, j, 0) - same_orient*valeur ;
536
537 valeur = sol_deux(2).val_point(air, theta, phi) ;
538 limite_z_un.set(num_front, k, j, 0) =
539 bound_z_un(num_front, k, j, 0) - valeur ;
540 }
541
542 // On modifie les Cl sur sol_deux :
543 for (int k=0 ; k<nbrep_deux ; k++)
544 for (int j=0 ; j<nbret_deux ; j++) {
545 xabs = xa_mtbl_deux (num_front+1, k, j, 0) ;
546 yabs = ya_mtbl_deux (num_front+1, k, j, 0) ;
547 zabs = za_mtbl_deux (num_front+1, k, j, 0) ;
548
549 source_un.get_mp()->convert_absolute
550 (xabs, yabs, zabs, air, theta, phi) ;
551
552 valeur = sol_un(0).val_point(air, theta, phi) ;
553 limite_x_deux.set(num_front, k, j, 0) =
554 bound_x_deux(num_front, k, j, 0) - same_orient*valeur ;
555
556 valeur = sol_un(1).val_point(air, theta, phi) ;
557 limite_y_deux.set(num_front, k, j, 0) =
558 bound_y_deux(num_front, k, j, 0) - same_orient*valeur ;
559
560 valeur = sol_un(2).val_point(air, theta, phi) ;
561 limite_z_deux.set(num_front, k, j, 0) =
562 bound_z_deux(num_front, k, j, 0) - valeur ;
563 }
564
565 double erreur = 0 ;
566
567 for (int i=0 ; i<3 ; i++) {
568 Tbl diff_un (diffrelmax (sol_un_old(i), sol_un(i))) ;
569 for (int j=num_front+1 ; j<nz_un ; j++)
570 if (erreur<diff_un(j))
571 erreur = diff_un(j) ;
572 }
573
574 for (int i=0 ; i<3 ; i++) {
575 Tbl diff_deux (diffrelmax (sol_deux_old(i), sol_deux(i))) ;
576 for (int j=num_front+1 ; j<nz_deux ; j++)
577 if (erreur<diff_deux(j))
578 erreur = diff_deux(j) ;
579 }
580
581 cout << "Pas " << conte << " : Difference " << erreur << endl ;
582
583 if (erreur < precision)
584 indic = -1 ;
585 conte ++ ;
586 }
587}
588void poisson_vect_binaire ( double lambda,
589 const Vector& source_un, const Vector& source_deux,
590 const Valeur& bound_x_un, const Valeur& bound_y_un,
591 const Valeur& bound_z_un, const Valeur& bound_x_deux,
592 const Valeur& bound_y_deux, const Valeur& bound_z_deux,
593 Vector& sol_un, Vector& sol_deux, int num_front, double precision) {
594
595 sol_un.change_triad(source_un.get_mp().get_bvect_cart()) ;
596 sol_deux.change_triad(source_deux.get_mp().get_bvect_cart()) ;
597
598 // Les bases des deux vecteurs doivent etre alignees ou non alignees :
599
600 assert (sol_un.get_mp() == source_un.get_mp()) ;
601 assert (sol_deux.get_mp() == source_deux.get_mp()) ;
602
603 double orientation_un = sol_un.get_mp().get_rot_phi() ;
604 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
605
606 double orientation_deux = sol_deux.get_mp().get_rot_phi() ;
607 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
608
609 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
610
611 Valeur limite_x_un (bound_x_un.get_mg()) ;
612 limite_x_un = bound_x_un ;
613 Valeur limite_y_un (bound_y_un.get_mg()) ;
614 limite_y_un = bound_y_un ;
615 Valeur limite_z_un (bound_z_un.get_mg()) ;
616 limite_z_un = bound_z_un ;
617
618 Valeur limite_x_deux (bound_x_deux.get_mg()) ;
619 limite_x_deux = bound_x_deux ;
620 Valeur limite_y_deux (bound_y_deux.get_mg()) ;
621 limite_y_deux = bound_y_deux ;
622 Valeur limite_z_deux (bound_z_deux.get_mg()) ;
623 limite_z_deux = bound_z_deux ;
624
625 Valeur limite_chi_un (bound_x_un.get_mg()) ;
626 limite_chi_un = 0 ;
627 limite_chi_un.std_base_scal() ;
628
629 Valeur limite_chi_deux (bound_x_deux.get_mg()) ;
630 limite_chi_deux = 0 ;
631 limite_chi_deux.std_base_scal() ;
632
633 Mtbl xa_mtbl_un (source_un.get_mp().get_mg()) ;
634 xa_mtbl_un.set_etat_qcq() ;
635 Mtbl ya_mtbl_un (source_un.get_mp().get_mg()) ;
636 ya_mtbl_un.set_etat_qcq() ;
637 Mtbl za_mtbl_un (source_un.get_mp().get_mg()) ;
638 za_mtbl_un.set_etat_qcq() ;
639 Mtbl xa_mtbl_deux (source_deux.get_mp().get_mg()) ;
640 xa_mtbl_deux.set_etat_qcq() ;
641 Mtbl ya_mtbl_deux (source_deux.get_mp().get_mg()) ;
642 ya_mtbl_deux.set_etat_qcq() ;
643 Mtbl za_mtbl_deux (source_deux.get_mp().get_mg()) ;
644 za_mtbl_deux.set_etat_qcq() ;
645
646 xa_mtbl_un = source_un.get_mp().xa ;
647 ya_mtbl_un = source_un.get_mp().ya ;
648 za_mtbl_un = source_un.get_mp().za ;
649
650 xa_mtbl_deux = source_deux.get_mp().xa ;
651 ya_mtbl_deux = source_deux.get_mp().ya ;
652 za_mtbl_deux = source_deux.get_mp().za ;
653
654 double xabs, yabs, zabs ;
655 double air, theta, phi ;
656 double valeur ;
657
658 int nbrep_un = bound_x_un.get_mg()->get_np(num_front) ;
659 int nbret_un = bound_x_un.get_mg()->get_nt(num_front) ;
660 int nbrep_deux = bound_x_deux.get_mg()->get_np(num_front) ;
661 int nbret_deux = bound_x_deux.get_mg()->get_nt(num_front) ;
662 int nz_un = bound_x_un.get_mg()->get_nzone() ;
663 int nz_deux = bound_x_deux.get_mg()->get_nzone() ;
664
665 const Metric_flat& ff_un (source_un.get_mp().flat_met_cart()) ;
666 const Metric_flat& ff_deux (source_deux.get_mp().flat_met_cart()) ;
667
668 // La source de l'equation scalaire sur 1
669 Vector cop_so_un (source_un) ;
670 cop_so_un.dec_dzpuis(2) ;
671 cop_so_un.dec_dzpuis(2) ;
672 cop_so_un.change_triad(source_un.get_mp().get_bvect_cart()) ;
673
674
675 Scalar source_scal_un (contract (cop_so_un.derive_cov(ff_un), 0, 1)/(lambda+1)) ;
676 if (source_scal_un.get_etat() == ETATZERO) {
677 source_scal_un.annule_hard() ;
678 source_scal_un.std_spectral_base() ;
679 }
680 source_scal_un.inc_dzpuis(2) ;
681
682 // La source de l'equation scalaire sur 2
683 Vector cop_so_deux (source_deux) ;
684 cop_so_deux.dec_dzpuis(2) ;
685 cop_so_deux.dec_dzpuis(2) ;
686 cop_so_deux.change_triad(source_deux.get_mp().get_bvect_cart()) ;
687
688 Scalar source_scal_deux (contract (cop_so_deux.derive_cov(ff_deux), 0, 1)/(lambda+1)) ;
689 if (source_scal_deux.get_etat() == ETATZERO) {
690 source_scal_deux.annule_hard() ;
691 source_scal_deux.std_spectral_base() ;
692 }
693 source_scal_deux.inc_dzpuis(2) ;
694
695 // Les copies :
696 Vector copie_so_un (source_un) ;
697 copie_so_un.dec_dzpuis() ;
698 copie_so_un.change_triad(source_un.get_mp().get_bvect_cart()) ;
699
700 Vector copie_so_deux (source_deux) ;
701 copie_so_deux.dec_dzpuis() ;
702 copie_so_deux.change_triad(source_deux.get_mp().get_bvect_cart()) ;
703
704 // ON COMMENCE LA BOUCLE :
705 Vector sol_un_old (sol_un) ;
706 Vector sol_deux_old (sol_deux) ;
707
708 int indic = 1 ;
709 int conte = 0 ;
710
711 while (indic == 1) {
712
713 // On resout les deux equations scalaires :
714 Scalar chi_un (source_scal_un.poisson_dirichlet (limite_chi_un, num_front)) ;
715 Scalar chi_deux (source_scal_deux.poisson_dirichlet (limite_chi_deux, num_front)) ;
716
717
718 // On calcule les sources pour les equations vectorielles :
719 Vector source_vect_un (copie_so_un) ;
720 Vector grad_chi_un (chi_un.derive_con(ff_un)) ;
721 grad_chi_un.inc_dzpuis() ;
722 source_vect_un = source_vect_un - lambda*grad_chi_un ;
723
724
725 Vector source_vect_deux (copie_so_deux) ;
726 Vector grad_chi_deux (chi_deux.derive_con(ff_deux)) ;
727 grad_chi_deux.inc_dzpuis() ;
728 source_vect_deux = source_vect_deux - lambda*grad_chi_deux ;
729
730 sol_un_old = sol_un ;
731 sol_deux_old = sol_deux ;
732
733 // On resout les equations vectorielles :
734 sol_un.set(1) = source_vect_un(1).poisson_dirichlet (limite_x_un, num_front) ;
735 sol_un.set(2) = source_vect_un(2).poisson_dirichlet (limite_y_un, num_front) ;
736 sol_un.set(3) = source_vect_un(3).poisson_dirichlet (limite_z_un, num_front) ;
737 sol_deux.set(1) = source_vect_deux(1).poisson_dirichlet (limite_x_deux, num_front) ;
738 sol_deux.set(2) = source_vect_deux(2).poisson_dirichlet (limite_y_deux, num_front) ;
739 sol_deux.set(3) = source_vect_deux(3).poisson_dirichlet (limite_z_deux, num_front) ;
740
741
742 // On modifie les Cl sur chi :
743 Scalar div_shift_un (contract(sol_un.derive_cov(ff_un), 0, 1)) ;
744 div_shift_un.dec_dzpuis(2) ;
745 div_shift_un.get_spectral_va().coef_i() ;
746
747 limite_chi_un = 1 ; // Affectation
748 for (int k=0 ; k<nbrep_un ; k++)
749 for (int j=0 ; j<nbret_un ; j++)
750 limite_chi_un.set(num_front, k, j, 0) =
751 div_shift_un.get_spectral_va() (num_front+1, k, j, 0) ;
752 limite_chi_un.std_base_scal() ;
753
754 Scalar div_shift_deux (contract(sol_deux.derive_cov(ff_deux), 0, 1)) ;
755 div_shift_deux.dec_dzpuis(2) ;
756 div_shift_deux.get_spectral_va().coef_i() ;
757
758 limite_chi_deux = 1 ; // Affectation
759 for (int k=0 ; k<nbrep_deux ; k++)
760 for (int j=0 ; j<nbret_deux ; j++)
761 limite_chi_deux.set(num_front, k, j, 0) =
762 div_shift_deux.get_spectral_va() (num_front+1, k, j, 0) ;
763 limite_chi_deux.std_base_scal() ;
764
765
766 // On modifie les Cl sur sol_un :
767 for (int k=0 ; k<nbrep_un ; k++)
768 for (int j=0 ; j<nbret_un ; j++) {
769 xabs = xa_mtbl_un (num_front+1, k, j, 0) ;
770 yabs = ya_mtbl_un (num_front+1, k, j, 0) ;
771 zabs = za_mtbl_un (num_front+1, k, j, 0) ;
772
773 source_deux.get_mp().convert_absolute
774 (xabs, yabs, zabs, air, theta, phi) ;
775
776 valeur = sol_deux(1).val_point(air, theta, phi) ;
777 limite_x_un.set(num_front, k, j, 0) =
778 bound_x_un(num_front, k, j, 0) - same_orient*valeur ;
779
780 valeur = sol_deux(2).val_point(air, theta, phi) ;
781 limite_y_un.set(num_front, k, j, 0) =
782 bound_y_un(num_front, k, j, 0) - same_orient*valeur ;
783
784 valeur = sol_deux(3).val_point(air, theta, phi) ;
785 limite_z_un.set(num_front, k, j, 0) =
786 bound_z_un(num_front, k, j, 0) - valeur ;
787 }
788
789 // On modifie les Cl sur sol_deux :
790 for (int k=0 ; k<nbrep_deux ; k++)
791 for (int j=0 ; j<nbret_deux ; j++) {
792 xabs = xa_mtbl_deux (num_front+1, k, j, 0) ;
793 yabs = ya_mtbl_deux (num_front+1, k, j, 0) ;
794 zabs = za_mtbl_deux (num_front+1, k, j, 0) ;
795
796 source_un.get_mp().convert_absolute
797 (xabs, yabs, zabs, air, theta, phi) ;
798
799 valeur = sol_un(1).val_point(air, theta, phi) ;
800 limite_x_deux.set(num_front, k, j, 0) =
801 bound_x_deux(num_front, k, j, 0) - same_orient*valeur ;
802
803 valeur = sol_un(2).val_point(air, theta, phi) ;
804 limite_y_deux.set(num_front, k, j, 0) =
805 bound_y_deux(num_front, k, j, 0) - same_orient*valeur ;
806
807 valeur = sol_un(3).val_point(air, theta, phi) ;
808 limite_z_deux.set(num_front, k, j, 0) =
809 bound_z_deux(num_front, k, j, 0) - valeur ;
810 }
811
812 double erreur = 0 ;
813
814 for (int i=1 ; i<=3 ; i++) {
815 Tbl diff_un (diffrelmax (sol_un_old(i), sol_un(i))) ;
816 for (int j=num_front+1 ; j<nz_un ; j++)
817 if (erreur<diff_un(j))
818 erreur = diff_un(j) ;
819 }
820
821 for (int i=1 ; i<=3 ; i++) {
822 Tbl diff_deux (diffrelmax (sol_deux_old(i), sol_deux(i))) ;
823 for (int j=num_front+1 ; j<nz_deux ; j++)
824 if (erreur<diff_deux(j))
825 erreur = diff_deux(j) ;
826 }
827
828 cout << "Pas " << conte << " : Difference " << erreur << endl ;
829
830
831
832/*
833 maxabs(sol_un.derive_con(ff_un).divergence(ff_un) + lambda * sol_un.divergence(ff_un).derive_con(ff_un) - source_un,
834 "Absolute error in the resolution of the equation for beta") ;
835 cout << endl ;
836*/
837
838
839 if (erreur < precision)
840 indic = -1 ;
841 conte ++ ;
842 }
843}
844}
void dec2_dzpuis()
dzpuis -= 2 ;
Definition tenseur.C:1130
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:481
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:539
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64