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 char 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
84 namespace Lorene {
85 void 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
197 void 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 
316 void 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 }
588 void 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 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:539
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:481
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Lorene::Tenseur::dec2_dzpuis
void dec2_dzpuis()
dzpuis -= 2 ;
Definition: tenseur.C:1130
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition: tenseur_operateur.C:279