LORENE
pde_frontiere_bin.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: pde_frontiere_bin.C,v 1.9 2018/11/16 14:34:36 j_novak Exp $
27  * $Log: pde_frontiere_bin.C,v $
28  * Revision 1.9 2018/11/16 14:34:36 j_novak
29  * Changed minor points to avoid some compilation warnings.
30  *
31  * Revision 1.8 2016/12/05 16:18:09 j_novak
32  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
33  *
34  * Revision 1.7 2014/10/13 08:53:29 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.6 2014/10/06 15:16:08 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.5 2006/05/11 14:16:37 f_limousin
41  * Minor modifs.
42  *
43  * Revision 1.4 2005/04/29 14:06:18 f_limousin
44  * Improve the resolution for Neumann_binaire(Scalars).
45  *
46  * Revision 1.3 2005/02/08 10:05:53 f_limousin
47  * Implementation of neumann_binaire(...) and dirichlet_binaire(...)
48  * with Scalars (instead of Cmps) in arguments.
49  *
50  * Revision 1.2 2003/10/03 15:58:50 j_novak
51  * Cleaning of some headers
52  *
53  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
54  * LORENE
55  *
56  * Revision 2.3 2000/12/04 14:30:16 phil
57  * correction CL.
58  *
59  * Revision 2.2 2000/12/01 15:18:26 phil
60  * vire trucs inutiles
61  *
62  * Revision 2.1 2000/12/01 15:16:49 phil
63  * correction version Neumann
64  *
65  * Revision 2.0 2000/10/19 09:36:07 phil
66  * *** empty log message ***
67  *
68  *
69  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/pde_frontiere_bin.C,v 1.9 2018/11/16 14:34:36 j_novak Exp $
70  *
71  */
72 
73 //standard
74 #include <cstdlib>
75 #include <cmath>
76 
77 // LORENE
78 #include "tensor.h"
79 #include "tenseur.h"
80 #include "proto.h"
81 #include "metric.h"
82 
83 // Version avec une fonction de theta, phi.
84 
85 namespace Lorene {
86 void dirichlet_binaire (const Cmp& source_un, const Cmp& source_deux,
87  const Valeur& boundary_un, const Valeur& boundary_deux,
88  Cmp& sol_un, Cmp& sol_deux, int num_front,
89  double precision) {
90 
91  // Les verifs sur le mapping :
92  assert (source_un.get_mp() == sol_un.get_mp()) ;
93  assert (source_deux.get_mp() == sol_deux.get_mp()) ;
94 
95  Valeur limite_un (boundary_un.get_mg()) ;
96  Valeur limite_deux (boundary_deux.get_mg()) ;
97 
98  Cmp sol_un_old (sol_un.get_mp()) ;
99  Cmp sol_deux_old (sol_deux.get_mp()) ;
100 
101  Mtbl xa_mtbl_un (source_un.get_mp()->xa) ;
102  Mtbl ya_mtbl_un (source_un.get_mp()->ya) ;
103  Mtbl za_mtbl_un (source_un.get_mp()->za) ;
104  Mtbl xa_mtbl_deux (source_deux.get_mp()->xa) ;
105  Mtbl ya_mtbl_deux (source_deux.get_mp()->ya) ;
106  Mtbl za_mtbl_deux (source_deux.get_mp()->za) ;
107 
108  double xabs, yabs, zabs ;
109  double air, theta, phi ;
110  double valeur ;
111 
112  int nbrep_un = boundary_un.get_mg()->get_np(num_front) ;
113  int nbret_un = boundary_un.get_mg()->get_nt(num_front) ;
114  int nbrep_deux = boundary_deux.get_mg()->get_np(num_front) ;
115  int nbret_deux = boundary_deux.get_mg()->get_nt(num_front) ;
116 
117  int nz_un = boundary_un.get_mg()->get_nzone() ;
118  int nz_deux = boundary_deux.get_mg()->get_nzone() ;
119 
120  // Initialisation valeur limite avant iteration !
121  limite_un = 1 ; //Pour initialiser les tableaux
122  for (int k=0 ; k<nbrep_un ; k++)
123  for (int j=0 ; j<nbret_un ; j++)
124  limite_un.set(num_front, k, j, 0) =
125  sol_un.va.val_point_jk(num_front+1, -1, j, k) ;
126  limite_un.set_base (boundary_un.base) ;
127 
128  limite_deux = 1 ;
129  for (int k=0 ; k<nbrep_deux ; k++)
130  for (int j=0 ; j<nbret_deux ; j++)
131  limite_deux.set(num_front, k, j, 0) =
132  sol_deux.va.val_point_jk(num_front+1, -1, j, k) ;
133  limite_deux.set_base (boundary_deux.base) ;
134 
135 
136  int conte = 0 ;
137  int indic = 1 ;
138 
139  while (indic==1) {
140 
141  sol_un_old = sol_un ;
142  sol_deux_old = sol_deux ;
143 
144  sol_un = source_un.poisson_dirichlet(limite_un, num_front) ;
145  sol_deux = source_deux.poisson_dirichlet(limite_deux, num_front) ;
146 
147  xa_mtbl_deux = source_deux.get_mp()->xa ;
148  ya_mtbl_deux = source_deux.get_mp()->ya ;
149  za_mtbl_deux = source_deux.get_mp()->za ;
150 
151 
152  for (int k=0 ; k<nbrep_deux ; k++)
153  for (int j=0 ; j<nbret_deux ; j++) {
154  xabs = xa_mtbl_deux (num_front+1, k, j, 0) ;
155  yabs = ya_mtbl_deux (num_front+1, k, j, 0) ;
156  zabs = za_mtbl_deux (num_front+1, k, j, 0) ;
157 
158  source_un.get_mp()->convert_absolute
159  (xabs, yabs, zabs, air, theta, phi) ;
160  valeur = sol_un.val_point(air, theta, phi) ;
161 
162  limite_deux.set(num_front, k, j, 0) =
163  boundary_deux(num_front, k, j, 0) - valeur ;
164  }
165 
166  xa_mtbl_un = source_un.get_mp()->xa ;
167  ya_mtbl_un = source_un.get_mp()->ya ;
168  za_mtbl_un = source_un.get_mp()->za ;
169 
170  for (int k=0 ; k<nbrep_un ; k++)
171  for (int j=0 ; j<nbret_un ; j++) {
172  xabs = xa_mtbl_un (num_front+1, k, j, 0) ;
173  yabs = ya_mtbl_un (num_front+1, k, j, 0) ;
174  zabs = za_mtbl_un (num_front+1, k, j, 0) ;
175 
176  source_deux.get_mp()->convert_absolute
177  (xabs, yabs, zabs, air, theta, phi) ;
178  valeur = sol_deux.val_point(air, theta, phi) ;
179 
180  limite_un.set(num_front, k, j, 0) =
181  boundary_un(num_front, k, j, 0) - valeur ;
182  }
183 
184  double erreur = 0 ;
185  Tbl diff_un (diffrelmax(sol_un, sol_un_old)) ;
186  for (int i=num_front+1 ; i<nz_un ; i++)
187  if (diff_un(i) > erreur)
188  erreur = diff_un(i) ;
189 
190  Tbl diff_deux (diffrelmax(sol_deux, sol_deux_old)) ;
191  for (int i=num_front+1 ; i<nz_deux ; i++)
192  if (diff_deux(i) > erreur)
193  erreur = diff_deux(i) ;
194 
195  cout << "Pas " << conte << " : Difference " << erreur << endl ;
196 
197  conte ++ ;
198  if (erreur < precision)
199  indic = -1 ;
200  }
201 
202 }
203 
204 // Version avec des doubles :
205 void dirichlet_binaire (const Cmp& source_un, const Cmp& source_deux,
206  double bound_un, double bound_deux,
207  Cmp& sol_un, Cmp& sol_deux, int num_front,
208  double precision) {
209 
210  Valeur boundary_un (source_un.get_mp()->get_mg()->get_angu()) ;
211  if (bound_un == 0)
212  boundary_un.annule_hard() ;
213  else
214  boundary_un = bound_un ;
215  boundary_un.std_base_scal() ;
216 
217  Valeur boundary_deux (source_deux.get_mp()->get_mg()->get_angu()) ;
218  if (bound_deux == 0)
219  boundary_deux.annule_hard() ;
220  else
221  boundary_deux = bound_deux ;
222  boundary_deux.std_base_scal() ;
223 
224  dirichlet_binaire (source_un, source_deux, boundary_un, boundary_deux,
225  sol_un, sol_deux, num_front, precision) ;
226 }
227 
228 
229 // Version with Scalar :
230 void dirichlet_binaire (const Scalar& source_un, const Scalar& source_deux,
231  const Valeur& boundary_un, const Valeur& boundary_deux,
232  Scalar& sol_un, Scalar& sol_deux, int num_front,
233  double precision) {
234 
235  // Les verifs sur le mapping :
236  assert (source_un.get_mp() == sol_un.get_mp()) ;
237  assert (source_deux.get_mp() == sol_deux.get_mp()) ;
238 
239  Valeur limite_un (boundary_un.get_mg()) ;
240  Valeur limite_deux (boundary_deux.get_mg()) ;
241 
242  Scalar sol_un_old (sol_un.get_mp()) ;
243  Scalar sol_deux_old (sol_deux.get_mp()) ;
244 
245  Mtbl xa_mtbl_un (source_un.get_mp().xa) ;
246  Mtbl ya_mtbl_un (source_un.get_mp().ya) ;
247  Mtbl za_mtbl_un (source_un.get_mp().za) ;
248  Mtbl xa_mtbl_deux (source_deux.get_mp().xa) ;
249  Mtbl ya_mtbl_deux (source_deux.get_mp().ya) ;
250  Mtbl za_mtbl_deux (source_deux.get_mp().za) ;
251 
252  double xabs, yabs, zabs ;
253  double air, theta, phi ;
254  double valeur ;
255 
256  int nbrep_un = boundary_un.get_mg()->get_np(num_front) ;
257  int nbret_un = boundary_un.get_mg()->get_nt(num_front) ;
258  int nbrep_deux = boundary_deux.get_mg()->get_np(num_front) ;
259  int nbret_deux = boundary_deux.get_mg()->get_nt(num_front) ;
260 
261  int nz_un = boundary_un.get_mg()->get_nzone() ;
262  int nz_deux = boundary_deux.get_mg()->get_nzone() ;
263 
264  // Initialisation valeur limite avant iteration !
265  limite_un = 1 ; //Pour initialiser les tableaux
266  for (int k=0 ; k<nbrep_un ; k++)
267  for (int j=0 ; j<nbret_un ; j++)
268  limite_un.set(num_front, k, j, 0) =
269  sol_un.get_spectral_va().val_point_jk(num_front+1, -1, j, k) ;
270  limite_un.set_base (boundary_un.base) ;
271 
272  limite_deux = 1 ;
273  for (int k=0 ; k<nbrep_deux ; k++)
274  for (int j=0 ; j<nbret_deux ; j++)
275  limite_deux.set(num_front, k, j, 0) =
276  sol_deux.get_spectral_va().val_point_jk(num_front+1, -1, j, k) ;
277  limite_deux.set_base (boundary_deux.base) ;
278 
279 
280  int conte = 0 ;
281  int indic = 1 ;
282 
283  while (indic==1) {
284 
285  sol_un_old = sol_un ;
286  sol_deux_old = sol_deux ;
287 
288  sol_un = source_un.poisson_dirichlet(limite_un, num_front) ;
289  sol_deux = source_deux.poisson_dirichlet(limite_deux, num_front) ;
290 
291  xa_mtbl_deux = source_deux.get_mp().xa ;
292  ya_mtbl_deux = source_deux.get_mp().ya ;
293  za_mtbl_deux = source_deux.get_mp().za ;
294 
295 
296  for (int k=0 ; k<nbrep_deux ; k++)
297  for (int j=0 ; j<nbret_deux ; j++) {
298  xabs = xa_mtbl_deux (num_front+1, k, j, 0) ;
299  yabs = ya_mtbl_deux (num_front+1, k, j, 0) ;
300  zabs = za_mtbl_deux (num_front+1, k, j, 0) ;
301 
302  source_un.get_mp().convert_absolute
303  (xabs, yabs, zabs, air, theta, phi) ;
304  valeur = sol_un.val_point(air, theta, phi) ;
305 
306  limite_deux.set(num_front, k, j, 0) =
307  boundary_deux(num_front, k, j, 0) - valeur ;
308  }
309 
310  xa_mtbl_un = source_un.get_mp().xa ;
311  ya_mtbl_un = source_un.get_mp().ya ;
312  za_mtbl_un = source_un.get_mp().za ;
313 
314  for (int k=0 ; k<nbrep_un ; k++)
315  for (int j=0 ; j<nbret_un ; j++) {
316  xabs = xa_mtbl_un (num_front+1, k, j, 0) ;
317  yabs = ya_mtbl_un (num_front+1, k, j, 0) ;
318  zabs = za_mtbl_un (num_front+1, k, j, 0) ;
319 
320  source_deux.get_mp().convert_absolute
321  (xabs, yabs, zabs, air, theta, phi) ;
322  valeur = sol_deux.val_point(air, theta, phi) ;
323 
324  limite_un.set(num_front, k, j, 0) =
325  boundary_un(num_front, k, j, 0) - valeur ;
326 // cout << 0.3 << " " << valeur+(sol_un+1.).val_grid_point(1, k, j, 0) << " " << (0.3 - (sol_un+1.)).val_grid_point(1, k, j, 0)-valeur << endl ;
327  }
328 
329  double erreur = 0 ;
330  Tbl diff_un (diffrelmax(sol_un, sol_un_old)) ;
331  for (int i=num_front+1 ; i<nz_un ; i++)
332  if (diff_un(i) > erreur)
333  erreur = diff_un(i) ;
334 
335  Tbl diff_deux (diffrelmax(sol_deux, sol_deux_old)) ;
336  for (int i=num_front+1 ; i<nz_deux ; i++)
337  if (diff_deux(i) > erreur)
338  erreur = diff_deux(i) ;
339 
340  cout << "Pas " << conte << " : Difference " << erreur << endl ;
341 
342  conte ++ ;
343  if (erreur < precision)
344  indic = -1 ;
345  }
346 
347 }
348 
349 
350 
351 // Version avec une fonction de theta, phi.
352 
353 void neumann_binaire (const Cmp& source_un, const Cmp& source_deux,
354  const Valeur& boundary_un, const Valeur& boundary_deux,
355  Cmp& sol_un, Cmp& sol_deux, int num_front,
356  double precision) {
357 
358  // Les verifs sur le mapping :
359  assert (source_un.get_mp() == sol_un.get_mp()) ;
360  assert (source_deux.get_mp() == sol_deux.get_mp()) ;
361 
362  // Alignes ou non ?
363  double orient_un = source_un.get_mp()->get_rot_phi() ;
364  assert ((orient_un==0) || (orient_un==M_PI)) ;
365  double orient_deux = source_deux.get_mp()->get_rot_phi() ;
366  assert ((orient_deux==0) || (orient_deux==M_PI)) ;
367  int same_orient = (orient_un==orient_deux) ? 1 : -1 ;
368 
369  Valeur limite_un (boundary_un.get_mg()) ;
370  Valeur limite_deux (boundary_deux.get_mg()) ;
371 
372  Cmp sol_un_old (sol_un.get_mp()) ;
373  Cmp sol_deux_old (sol_deux.get_mp()) ;
374 
375  Mtbl xa_mtbl_un (source_un.get_mp()->xa) ;
376  Mtbl ya_mtbl_un (source_un.get_mp()->ya) ;
377  Mtbl za_mtbl_un (source_un.get_mp()->za) ;
378 
379  Mtbl cost_mtbl_un (source_un.get_mp()->cost) ;
380  Mtbl sint_mtbl_un (source_un.get_mp()->sint) ;
381  Mtbl cosp_mtbl_un (source_un.get_mp()->cosp) ;
382  Mtbl sinp_mtbl_un (source_un.get_mp()->sinp) ;
383 
384 
385  Mtbl xa_mtbl_deux (source_deux.get_mp()->xa) ;
386  Mtbl ya_mtbl_deux (source_deux.get_mp()->ya) ;
387  Mtbl za_mtbl_deux (source_deux.get_mp()->za) ;
388 
389  Mtbl cost_mtbl_deux (source_deux.get_mp()->cost) ;
390  Mtbl sint_mtbl_deux (source_deux.get_mp()->sint) ;
391  Mtbl cosp_mtbl_deux (source_deux.get_mp()->cosp) ;
392  Mtbl sinp_mtbl_deux (source_deux.get_mp()->sinp) ;
393 
394  double xabs, yabs, zabs ;
395  double air, theta, phi ;
396  double valeur ;
397 
398  int nbrep_un = boundary_un.get_mg()->get_np(num_front) ;
399  int nbret_un = boundary_un.get_mg()->get_nt(num_front) ;
400  int nbrep_deux = boundary_deux.get_mg()->get_np(num_front) ;
401  int nbret_deux = boundary_deux.get_mg()->get_nt(num_front) ;
402 
403  int nz_un = boundary_un.get_mg()->get_nzone() ;
404  int nz_deux = boundary_deux.get_mg()->get_nzone() ;
405 
406  // Initialisation des CL :
407  limite_un = 1 ;
408  limite_deux = 2 ;
409  Cmp der_un (sol_un.dsdr()) ;
410  Cmp der_deux (sol_deux.dsdr()) ;
411 
412  for (int k=0 ; k<nbrep_un ; k++)
413  for (int j=0 ; j<nbret_un ; j++)
414  limite_un.set(num_front, k, j, 0) =
415  der_un.va.val_point_jk(num_front+1, -1, j, k) ;
416  limite_un.set_base (boundary_un.base) ;
417 
418  for (int k=0 ; k<nbrep_deux ; k++)
419  for (int j=0 ; j<nbret_deux ; j++)
420  limite_deux.set(num_front, k, j, 0) =
421  der_deux.va.val_point_jk(num_front+1, -1, j, k) ;
422  limite_deux.set_base (boundary_deux.base) ;
423 
424  int conte = 0 ;
425  int indic = 1 ;
426 
427  while (indic==1) {
428 
429  sol_un_old = sol_un ;
430  sol_deux_old = sol_deux ;
431 
432  sol_un = source_un.poisson_neumann(limite_un, num_front) ;
433  sol_deux = source_deux.poisson_neumann(limite_deux, num_front) ;
434 
435  // On veut les derivees de l'un sur l'autre :
436  Tenseur copie_un (sol_un) ;
437  Tenseur grad_sol_un (copie_un.gradient()) ;
438  grad_sol_un.dec2_dzpuis() ;
439  grad_sol_un.set(0) = grad_sol_un(0)*same_orient ;
440  grad_sol_un.set(1) = grad_sol_un(1)*same_orient ;
441 
442  for (int k=0 ; k<nbrep_deux ; k++)
443  for (int j=0 ; j<nbret_deux ; j++) {
444  xabs = xa_mtbl_deux (num_front+1, k, j, 0) ;
445  yabs = ya_mtbl_deux (num_front+1, k, j, 0) ;
446  zabs = za_mtbl_deux (num_front+1, k, j, 0) ;
447 
448  source_un.get_mp()->convert_absolute
449  (xabs, yabs, zabs, air, theta, phi) ;
450 
451  valeur = sint_mtbl_deux (num_front+1, k, j, 0) * (
452 cosp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(0).val_point(air, theta, phi)+
453 sinp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(1).val_point(air, theta, phi))+
454 cost_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(2).val_point(air, theta, phi);
455 
456  limite_deux.set(num_front, k, j, 0) =
457  boundary_deux(num_front, k, j, 0) - valeur ;
458  }
459 
460  Tenseur copie_deux (sol_deux) ;
461  Tenseur grad_sol_deux (copie_deux.gradient()) ;
462  grad_sol_deux.dec2_dzpuis() ;
463  grad_sol_deux.set(0) = grad_sol_deux(0)*same_orient ;
464  grad_sol_deux.set(1) = grad_sol_deux(1)*same_orient ;
465 
466  for (int k=0 ; k<nbrep_un ; k++)
467  for (int j=0 ; j<nbret_un ; j++) {
468  xabs = xa_mtbl_un (num_front+1, k, j, 0) ;
469  yabs = ya_mtbl_un (num_front+1, k, j, 0) ;
470  zabs = za_mtbl_un (num_front+1, k, j, 0) ;
471 
472  source_deux.get_mp()->convert_absolute
473  (xabs, yabs, zabs, air, theta, phi) ;
474 
475  valeur = sint_mtbl_un (num_front+1, k, j, 0) * (
476 cosp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(0).val_point(air, theta, phi)+
477 sinp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(1).val_point(air, theta, phi))+
478 cost_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(2).val_point(air, theta, phi);
479 
480  limite_un.set(num_front, k, j, 0) =
481  boundary_un(num_front, k, j, 0) - valeur ;
482  }
483 
484  double erreur = 0 ;
485  Tbl diff_un (diffrelmax(sol_un, sol_un_old)) ;
486  for (int i=num_front+1 ; i<nz_un ; i++)
487  if (diff_un(i) > erreur)
488  erreur = diff_un(i) ;
489 
490  Tbl diff_deux (diffrelmax(sol_deux, sol_deux_old)) ;
491  for (int i=num_front+1 ; i<nz_deux ; i++)
492  if (diff_deux(i) > erreur)
493  erreur = diff_deux(i) ;
494 
495  cout << "Pas " << conte << " : Difference " << erreur << endl ;
496  conte ++ ;
497 
498  if (erreur < precision)
499  indic = -1 ;
500  }
501 }
502 
503 // Version avec des doubles :
504 void neumann_binaire (const Cmp& source_un, const Cmp& source_deux,
505  double bound_un, double bound_deux,
506  Cmp& sol_un, Cmp& sol_deux, int num_front,
507  double precision) {
508 
509  Valeur boundary_un (source_un.get_mp()->get_mg()->get_angu()) ;
510  if (bound_un == 0)
511  boundary_un.annule_hard () ;
512  else
513  boundary_un = bound_un ;
514  boundary_un.std_base_scal() ;
515 
516  Valeur boundary_deux (source_deux.get_mp()->get_mg()->get_angu()) ;
517  if (bound_deux == 0)
518  boundary_deux.annule_hard() ;
519  else
520  boundary_deux = bound_deux ;
521  boundary_deux.std_base_scal() ;
522 
523  neumann_binaire (source_un, source_deux, boundary_un, boundary_deux,
524  sol_un, sol_deux, num_front, precision) ;
525 }
526 void neumann_binaire (const Scalar& source_un, const Scalar& source_deux,
527  const Valeur& boundary_un, const Valeur& boundary_deux,
528  Scalar& sol_un, Scalar& sol_deux, int num_front,
529  double precision) {
530 
531  // Les verifs sur le mapping :
532  assert (source_un.get_mp() == sol_un.get_mp()) ;
533  assert (source_deux.get_mp() == sol_deux.get_mp()) ;
534 
535  // Alignement
536  double orient_un = source_un.get_mp().get_rot_phi() ;
537  assert ((orient_un==0) || (orient_un==M_PI)) ;
538  double orient_deux = source_deux.get_mp().get_rot_phi() ;
539  assert ((orient_deux==0) || (orient_deux==M_PI)) ;
540  int same_orient = (orient_un==orient_deux) ? 1 : -1 ;
541 
542  Valeur limite_un (boundary_un.get_mg()) ;
543  Valeur limite_deux (boundary_deux.get_mg()) ;
544 
545  Scalar sol_un_old (sol_un.get_mp()) ;
546  Scalar sol_deux_old (sol_deux.get_mp()) ;
547 
548  Mtbl xa_mtbl_un (source_un.get_mp().xa) ;
549  Mtbl ya_mtbl_un (source_un.get_mp().ya) ;
550  Mtbl za_mtbl_un (source_un.get_mp().za) ;
551 
552  Mtbl cost_mtbl_un (source_un.get_mp().cost) ;
553  Mtbl sint_mtbl_un (source_un.get_mp().sint) ;
554  Mtbl cosp_mtbl_un (source_un.get_mp().cosp) ;
555  Mtbl sinp_mtbl_un (source_un.get_mp().sinp) ;
556 
557  Mtbl xa_mtbl_deux (source_deux.get_mp().xa) ;
558  Mtbl ya_mtbl_deux (source_deux.get_mp().ya) ;
559  Mtbl za_mtbl_deux (source_deux.get_mp().za) ;
560 
561  Mtbl cost_mtbl_deux (source_deux.get_mp().cost) ;
562  Mtbl sint_mtbl_deux (source_deux.get_mp().sint) ;
563  Mtbl cosp_mtbl_deux (source_deux.get_mp().cosp) ;
564  Mtbl sinp_mtbl_deux (source_deux.get_mp().sinp) ;
565 
566  double xabs, yabs, zabs ;
567  double air, theta, phi ;
568  double valeur ;
569 
570  const Metric_flat& ff_un (source_un.get_mp().flat_met_cart()) ;
571  const Metric_flat& ff_deux (source_deux.get_mp().flat_met_cart()) ;
572 
573  int nbrep_un = boundary_un.get_mg()->get_np(num_front) ;
574  int nbret_un = boundary_un.get_mg()->get_nt(num_front) ;
575  int nbrep_deux = boundary_deux.get_mg()->get_np(num_front) ;
576  int nbret_deux = boundary_deux.get_mg()->get_nt(num_front) ;
577 
578  int nz_un = boundary_un.get_mg()->get_nzone() ;
579  int nz_deux = boundary_deux.get_mg()->get_nzone() ;
580 
581  // Initialisation des CL :
582  limite_un = 1 ;
583  limite_deux = 2 ;
584  Scalar der_un (sol_un.dsdr()) ;
585  Scalar der_deux (sol_deux.dsdr()) ;
586 
587  for (int k=0 ; k<nbrep_un ; k++)
588  for (int j=0 ; j<nbret_un ; j++)
589  limite_un.set(num_front, k, j, 0) =
590  der_un.get_spectral_va().val_point_jk(num_front+1, -1, j, k) ;
591  limite_un.set_base (boundary_un.base) ;
592 
593  for (int k=0 ; k<nbrep_deux ; k++)
594  for (int j=0 ; j<nbret_deux ; j++)
595  limite_deux.set(num_front, k, j, 0) =
596  der_deux.get_spectral_va().val_point_jk(num_front+1, -1, j, k) ;
597  limite_deux.set_base (boundary_deux.base) ;
598 
599  int conte = 0 ;
600  int indic = 1 ;
601 
602  while (indic==1) {
603 
604  sol_un_old = sol_un ;
605  sol_deux_old = sol_deux ;
606 
607  sol_un = source_un.poisson_neumann(limite_un, num_front) ;
608  sol_deux = source_deux.poisson_neumann(limite_deux, num_front) ;
609 
610  // On veut les derivees de l'un sur l'autre :
611  Scalar copie_un (sol_un) ;
612  Vector grad_sol_un (copie_un.derive_cov(ff_un)) ;
613  grad_sol_un.dec_dzpuis(2) ;
614  grad_sol_un.set(1) = grad_sol_un(1)*same_orient ;
615  grad_sol_un.set(2) = grad_sol_un(2)*same_orient ;
616 
617 
618  for (int k=0 ; k<nbrep_deux ; k++)
619  for (int j=0 ; j<nbret_deux ; j++) {
620  xabs = xa_mtbl_deux (num_front+1, k, j, 0) ;
621  yabs = ya_mtbl_deux (num_front+1, k, j, 0) ;
622  zabs = za_mtbl_deux (num_front+1, k, j, 0) ;
623 
624  source_un.get_mp().convert_absolute
625  (xabs, yabs, zabs, air, theta, phi) ;
626 
627  valeur = sint_mtbl_deux (num_front+1, k, j, 0) * (
628 cosp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(1).val_point(air, theta, phi)+
629 sinp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(2).val_point(air, theta, phi))+
630 cost_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(3).val_point(air, theta, phi);
631 
632  limite_deux.set(num_front, k, j, 0) =
633  boundary_deux(num_front, k, j, 0) - valeur ;
634  }
635 
636  Scalar copie_deux (sol_deux) ;
637  Vector grad_sol_deux (copie_deux.derive_cov(ff_deux)) ;
638  grad_sol_deux.dec_dzpuis(2) ;
639  grad_sol_deux.set(1) = grad_sol_deux(1)*same_orient ;
640  grad_sol_deux.set(2) = grad_sol_deux(2)*same_orient ;
641 
642  for (int k=0 ; k<nbrep_un ; k++)
643  for (int j=0 ; j<nbret_un ; j++) {
644  xabs = xa_mtbl_un (num_front+1, k, j, 0) ;
645  yabs = ya_mtbl_un (num_front+1, k, j, 0) ;
646  zabs = za_mtbl_un (num_front+1, k, j, 0) ;
647 
648  source_deux.get_mp().convert_absolute
649  (xabs, yabs, zabs, air, theta, phi) ;
650 
651  valeur = sint_mtbl_un (num_front+1, k, j, 0) * (
652 cosp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(1).val_point(air, theta, phi)+
653 sinp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(2).val_point(air, theta, phi))+
654 cost_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(3).val_point(air, theta, phi);
655 
656 
657 
658  limite_un.set(num_front, k, j, 0) =
659  boundary_un(num_front, k, j, 0) - valeur ;
660  }
661 
662  double erreur = 0 ;
663  Tbl diff_un (diffrelmax(sol_un, sol_un_old)) ;
664  for (int i=num_front+1 ; i<nz_un ; i++)
665  if (diff_un(i) > erreur)
666  erreur = diff_un(i) ;
667 
668  Tbl diff_deux (diffrelmax(sol_deux, sol_deux_old)) ;
669  for (int i=num_front+1 ; i<nz_deux ; i++)
670  if (diff_deux(i) > erreur)
671  erreur = diff_deux(i) ;
672 
673  cout << "Pas " << conte << " : Difference " << erreur << endl ;
674  conte ++ ;
675 
676  Scalar source1 (source_un) ;
677  Scalar solution1 (sol_un) ;
678 
679 // maxabs(solution1.laplacian() - source1,
680 // "Absolute error in the resolution of the equation for psi") ;
681 
682  if (erreur < precision)
683  indic = -1 ;
684  }
685 }
686 }
Lorene
Lorene prototypes.
Definition: app_hor.h:67
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:542