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