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
85namespace Lorene {
86void 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 :
205void 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 :
230void 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
353void 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) * (
452cosp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(0).val_point(air, theta, phi)+
453sinp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(1).val_point(air, theta, phi))+
454cost_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) * (
476cosp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(0).val_point(air, theta, phi)+
477sinp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(1).val_point(air, theta, phi))+
478cost_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 :
504void 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}
526void 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) * (
628cosp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(1).val_point(air, theta, phi)+
629sinp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(2).val_point(air, theta, phi))+
630cost_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) * (
652cosp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(1).val_point(air, theta, phi)+
653sinp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(2).val_point(air, theta, phi))+
654cost_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::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::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::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:542
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:732