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
23char 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
79namespace Lorene {
80void 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 :
199void 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 :
224void 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
347void 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) * (
446cosp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(0).val_point(air, theta, phi)+
447sinp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(1).val_point(air, theta, phi))+
448cost_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) * (
470cosp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(0).val_point(air, theta, phi)+
471sinp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(1).val_point(air, theta, phi))+
472cost_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 :
498void 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}
520void 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) * (
622cosp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(1).val_point(air, theta, phi)+
623sinp_mtbl_deux(num_front+1, k, j, 0)*grad_sol_un(2).val_point(air, theta, phi))+
624cost_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) * (
646cosp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(1).val_point(air, theta, phi)+
647sinp_mtbl_un(num_front+1, k, j, 0)*grad_sol_deux(2).val_point(air, theta, phi))+
648cost_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}
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:539
Lorene prototypes.
Definition app_hor.h:64