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