LORENE
sol_elliptic_sin_zec.C
1 /*
2  * Copyright (c) 2004 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 version 2
8  * as published by the Free Software Foundation.
9  *
10  * LORENE is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with LORENE; if not, write to the Free Software
17  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18  *
19  */
20 /*
21  * $Id: sol_elliptic_sin_zec.C,v 1.8 2014/10/13 08:53:30 j_novak Exp $
22  * $Log: sol_elliptic_sin_zec.C,v $
23  * Revision 1.8 2014/10/13 08:53:30 j_novak
24  * Lorene classes and functions now belong to the namespace Lorene.
25  *
26  * Revision 1.7 2014/10/06 15:16:10 j_novak
27  * Modified #include directives to use c++ syntax.
28  *
29  * Revision 1.6 2007/05/08 07:08:30 p_grandclement
30  * *** empty log message ***
31  *
32  * Revision 1.5 2007/05/06 10:48:12 p_grandclement
33  * Modification of a few operators for the vorton project
34  *
35  * Revision 1.4 2005/11/30 11:09:08 p_grandclement
36  * Changes for the Bin_ns_bh project
37  *
38  * Revision 1.3 2005/08/26 14:02:41 p_grandclement
39  * Modification of the elliptic solver that matches with an oscillatory exterior solution
40  * small correction in Poisson tau also...
41  *
42  * Revision 1.2 2004/08/24 09:14:44 p_grandclement
43  * Addition of some new operators, like Poisson in 2d... It now requieres the
44  * GSL library to work.
45  *
46  * Also, the way a variable change is stored by a Param_elliptic is changed and
47  * no longer uses Change_var but rather 2 Scalars. The codes using that feature
48  * will requiere some modification. (It should concern only the ones about monopoles)
49  *
50  * Revision 1.1 2004/02/11 09:52:52 p_grandclement
51  * Forgot one new file ...
52  *
53 
54  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/sol_elliptic_sin_zec.C,v 1.8 2014/10/13 08:53:30 j_novak Exp $
55  */
56 
57 char sol_elliptic_sin_zec_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/sol_elliptic_sin_zec.C,v 1.8 2014/10/13 08:53:30 j_novak Exp $" ;
58 
59 // Header C :
60 #include <cstdlib>
61 #include <cmath>
62 #include <gsl/gsl_sf_bessel.h>
63 
64 // Headers Lorene :
65 #include "tbl.h"
66 #include "mtbl_cf.h"
67 #include "map.h"
68 #include "param_elliptic.h"
69 
70 
71  //----------------------------------------------
72  // Version Mtbl_cf
73  //----------------------------------------------
74 
75 namespace Lorene {
76 Mtbl_cf elliptic_solver_sin_zec (const Param_elliptic& ope_var,
77  const Mtbl_cf& source, double* amplis, double* phases) {
78 
79  // Verifications d'usage sur les zones
80  int nz = source.get_mg()->get_nzone() ;
81  assert (nz>1) ;
82  assert (source.get_mg()->get_type_r(0) == RARE) ;
83  assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
84  for (int l=1 ; l<nz-1 ; l++)
85  assert(source.get_mg()->get_type_r(l) == FIN) ;
86 
87  // donnees sur la zone
88  int nr, nt, np ;
89 
90  //Rangement des valeurs intermediaires
91  Tbl *so ;
92  Tbl *sol_hom ;
93  Tbl *sol_part ;
94 
95 
96  // Rangement des solutions, avant raccordement
97  Mtbl_cf solution_part(source.get_mg(), source.base) ;
98  Mtbl_cf solution_hom_un(source.get_mg(), source.base) ;
99  Mtbl_cf solution_hom_deux(source.get_mg(), source.base) ;
100  Mtbl_cf resultat(source.get_mg(), source.base) ;
101 
102  solution_part.annule_hard() ;
103  solution_hom_un.annule_hard() ;
104  solution_hom_deux.annule_hard() ;
105  resultat.annule_hard() ;
106 
107  // Computation of the SP and SH's in every domain but the ZEC ...
108  int conte = 0 ;
109  for (int zone=0 ; zone<nz-1 ; zone++) {
110  nr = source.get_mg()->get_nr(zone) ;
111  nt = source.get_mg()->get_nt(zone) ;
112  np = source.get_mg()->get_np(zone) ;
113 
114  for (int k=0 ; k<np+1 ; k++)
115  for (int j=0 ; j<nt ; j++) {
116  if (ope_var.operateurs[conte] != 0x0) {
117  // Calcul de la SH
118  sol_hom = new Tbl(ope_var.operateurs[conte]->get_solh()) ;
119 
120  //Calcul de la SP
121  so = new Tbl(nr) ;
122  so->set_etat_qcq() ;
123  for (int i=0 ; i<nr ; i++)
124  so->set(i) = source(zone, k, j, i) ;
125 
126  sol_part = new Tbl(ope_var.operateurs[conte]->get_solp(*so)) ;
127 
128  // Rangement dans les tableaux globaux ;
129  for (int i=0 ; i<nr ; i++) {
130  solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
131  if (sol_hom->get_ndim()==1)
132  solution_hom_un.set(zone, k, j, i) = (*sol_hom)(i) ;
133  else
134  {
135  solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0,i) ;
136  solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1,i) ;
137  }
138  }
139 
140  delete so ;
141  delete sol_hom ;
142  delete sol_part ;
143 
144  }
145  conte ++ ;
146  }
147  }
148 
149  //-------------------------------------------------
150  // ON EST PARTI POUR LE RACCORD (Be carefull ....)
151  //-------------------------------------------------
152 
153  // C'est pas simple toute cette sombre affaire...
154  // Que le cas meme nombre de points dans chaque domaines...
155 
156  int start = 0 ;
157  for (int k=0 ; k< source.get_mg()->get_np(0)+1; k++)
158  for (int j=0 ; j<source.get_mg()->get_nt(0) ; j++) {
159  if (ope_var.operateurs[start] != 0x0) {
160 
161  int taille = 2*nz - 2 ;
162  Matrice systeme (taille, taille) ;
163  systeme.set_etat_qcq() ;
164  for (int i=0 ; i<taille ; i++)
165  for (int j2=0 ; j2<taille ; j2++)
166  systeme.set(i,j2) = 0 ;
167  Tbl sec_membre (taille) ;
168  sec_membre.set_etat_qcq() ;
169  for (int i=0 ; i<taille ; i++)
170  sec_membre.set(i) = 0 ;
171  //---------
172  // Noyau :
173  //---------
174  conte = start ;
175 
176  systeme.set(0,0) = ope_var.G_plus(0) *
177  ope_var.operateurs[conte]->val_sh_one_plus() ;
178  systeme.set(1,0) =
179  ope_var.dG_plus(0) * ope_var.operateurs[conte]->val_sh_one_plus() +
180  ope_var.G_plus(0) * ope_var.operateurs[conte]->der_sh_one_plus() ;
181 
182  sec_membre.set(0) -= ope_var.F_plus(0,k,j) +
183  ope_var.G_plus(0) * ope_var.operateurs[conte]->val_sp_plus() ;
184  sec_membre.set(1) -= ope_var.dF_plus(0,k,j) +
185  ope_var.dG_plus(0) * ope_var.operateurs[conte]->val_sp_plus() +
186  ope_var.G_plus(0) * ope_var.operateurs[conte]->der_sp_plus() ;
187 
188  //----------
189  // SHELLS :
190  //----------
191 
192  for (int l=1 ; l<nz-1 ; l++) {
193 
194  // On se met au bon endroit :
195  int np_prec = source.get_mg()->get_np(l-1) ;
196  int nt_prec = source.get_mg()->get_nt(l-1) ;
197  conte += (np_prec+1)*nt_prec ;
198 
199  systeme.set(2*l-2, 2*l-1) = -ope_var.G_minus(l) *
200  ope_var.operateurs[conte]->val_sh_one_minus() ;
201  systeme.set(2*l-2, 2*l) = - ope_var.G_minus(l) *
202  ope_var.operateurs[conte]->val_sh_two_minus() ;
203  systeme.set(2*l-1, 2*l-1) =
204  -ope_var.dG_minus(l)*ope_var.operateurs[conte]->val_sh_one_minus()-
205  ope_var.G_minus(l)*ope_var.operateurs[conte]->der_sh_one_minus() ;
206  systeme.set(2*l-1, 2*l) =
207  -ope_var.dG_minus(l)*ope_var.operateurs[conte]->val_sh_two_minus()-
208  ope_var.G_minus(l)*ope_var.operateurs[conte]->der_sh_two_minus() ;
209 
210  sec_membre.set(2*l-2) += ope_var.F_minus(l,k,j) +
211  ope_var.G_minus(l) * ope_var.operateurs[conte]->val_sp_minus() ;
212  sec_membre.set(2*l-1) += ope_var.dF_minus(l,k,j) +
213  ope_var.dG_minus(l) * ope_var.operateurs[conte]->val_sp_minus() +
214  ope_var.G_minus(l) * ope_var.operateurs[conte]->der_sp_minus() ;
215 
216  // Valeurs en +1 :
217  systeme.set(2*l, 2*l-1) = ope_var.G_plus(l) *
218  ope_var.operateurs[conte]->val_sh_one_plus() ;
219  systeme.set(2*l, 2*l) = ope_var.G_plus(l) *
220  ope_var.operateurs[conte]->val_sh_two_plus() ;
221 
222  systeme.set(2*l+1, 2*l-1) =
223  ope_var.dG_plus(l)*ope_var.operateurs[conte]->val_sh_one_plus()+
224  ope_var.G_plus(l)*ope_var.operateurs[conte]->der_sh_one_plus() ;
225  systeme.set(2*l+1, 2*l) =
226  ope_var.dG_plus(l)*ope_var.operateurs[conte]->val_sh_two_plus()+
227  ope_var.G_plus(l)*ope_var.operateurs[conte]->der_sh_two_plus() ;
228 
229  sec_membre.set(2*l) -= ope_var.F_plus(l,k,j) +
230  ope_var.G_plus(l) * ope_var.operateurs[conte]->val_sp_plus();
231 
232  sec_membre.set(2*l+1) -= ope_var.dF_plus(l,k,j) +
233  ope_var.dG_plus(l) * ope_var.operateurs[conte]->val_sp_plus() +
234  ope_var.G_plus(l) * ope_var.operateurs[conte]->der_sp_plus() ;
235 
236  }
237 
238  // On recupere la valeur de la sh :
239  double val_sh = cos(phases[start])*ope_var.operateurs[conte]->val_sh_one_plus()
240  + sin(phases[start])*ope_var.operateurs[conte]->val_sh_two_plus() ;
241  double der_sh = cos(phases[start])*ope_var.operateurs[conte]->der_sh_one_plus()
242  + sin(phases[start])*ope_var.operateurs[conte]->der_sh_two_plus() ;
243 
244  systeme.set(taille-2, taille-1) = -ope_var.G_minus(nz-1) * val_sh ;
245  systeme.set(taille-1, taille-1) = -ope_var.dG_minus(nz-1)*val_sh- ope_var.G_minus(nz-1)*der_sh ;
246 
247  sec_membre.set(taille-2) += ope_var.F_minus(nz-1,k,j) ;
248  sec_membre.set(taille-1) += ope_var.dF_minus(nz-1,k,j) ;
249 
250  // On resout le systeme ...
251  if (taille > 2)
252  systeme.set_band(2,2) ;
253  else
254  systeme.set_band(1,1) ;
255 
256  systeme.set_lu() ;
257  Tbl facteur (systeme.inverse(sec_membre)) ;
258 
259  amplis[start] = facteur(taille-1) ;
260 
261  // On range tout ca :
262  // Noyau
263  nr = source.get_mg()->get_nr(0) ;
264  for (int i=0 ; i<nr ; i++)
265  resultat.set(0,k,j,i) = solution_part(0,k,j,i)
266  +facteur(0)*solution_hom_un(0,k,j,i) ;
267 
268  // Shells
269  for (int l=1 ; l<nz-1 ; l++) {
270  nr = source.get_mg()->get_nr(l) ;
271  for (int i=0 ; i<nr ; i++)
272  resultat.set(l,k,j,i) = solution_part(l,k,j,i) +
273  facteur(2*l-1)*solution_hom_un(l,k,j,i) +
274  facteur(2*l)*solution_hom_deux(l,k,j,i) ;
275  }
276  }
277  start ++ ;
278  }
279  return resultat;
280 }
281 
282 
283 }
Lorene prototypes.
Definition: app_hor.h:64
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition: mtbl_cf.h:453
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69