LORENE
sol_elliptic.C
1 /*
2  * Copyright (c) 2003 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 char sol_elliptic_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/sol_elliptic.C,v 1.10 2014/10/13 08:53:30 j_novak Exp $" ;
22 
23 /*
24  * $Id: sol_elliptic.C,v 1.10 2014/10/13 08:53:30 j_novak Exp $
25  * $Log: sol_elliptic.C,v $
26  * Revision 1.10 2014/10/13 08:53:30 j_novak
27  * Lorene classes and functions now belong to the namespace Lorene.
28  *
29  * Revision 1.9 2014/10/06 15:16:10 j_novak
30  * Modified #include directives to use c++ syntax.
31  *
32  * Revision 1.8 2008/07/10 11:20:33 p_grandclement
33  * mistake fixed in solh_helmholtz_minus
34  *
35  * Revision 1.7 2008/07/09 06:51:58 p_grandclement
36  * some corrections to helmholtz minus in the nucleus
37  *
38  * Revision 1.6 2005/01/25 15:44:20 j_novak
39  * Fixed a problem in the definition of nr at the end.
40  *
41  * Revision 1.5 2004/11/25 08:14:56 j_novak
42  * Modifs. comments
43  *
44  * Revision 1.4 2004/08/24 09:14:44 p_grandclement
45  * Addition of some new operators, like Poisson in 2d... It now requieres the
46  * GSL library to work.
47  *
48  * Also, the way a variable change is stored by a Param_elliptic is changed and
49  * no longer uses Change_var but rather 2 Scalars. The codes using that feature
50  * will requiere some modification. (It should concern only the ones about monopoles)
51  *
52  * Revision 1.3 2004/05/26 20:32:44 p_grandclement
53  * utilisation of val_r_jk
54  *
55  * Revision 1.2 2003/12/19 16:21:49 j_novak
56  * Shadow hunt
57  *
58  * Revision 1.1 2003/12/11 14:48:49 p_grandclement
59  * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
60  *
61  *
62  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/sol_elliptic.C,v 1.10 2014/10/13 08:53:30 j_novak Exp $
63  *
64  */
65 
66 // Header C :
67 #include <cstdlib>
68 #include <cmath>
69 
70 // Headers Lorene :
71 #include "param_elliptic.h"
72 #include "tbl.h"
73 #include "mtbl_cf.h"
74 #include "map.h"
75 
76 
77  //----------------------------------------------
78  // Version Mtbl_cf
79  //----------------------------------------------
80 
81 
82 
83 namespace Lorene {
84 Mtbl_cf elliptic_solver (const Param_elliptic& ope_var, const Mtbl_cf& source) {
85  // Verifications d'usage sur les zones
86  int nz = source.get_mg()->get_nzone() ;
87  assert (nz>1) ;
88  assert (source.get_mg()->get_type_r(0) == RARE) ;
89  assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
90  for (int l=1 ; l<nz-1 ; l++)
91  assert(source.get_mg()->get_type_r(l) == FIN) ;
92 
93  // donnees sur la zone
94  int nr, nt, np ;
95 
96  //Rangement des valeurs intermediaires
97  Tbl *so ;
98  Tbl *sol_hom ;
99  Tbl *sol_part ;
100 
101 
102  // Rangement des solutions, avant raccordement
103  Mtbl_cf solution_part(source.get_mg(), source.base) ;
104  Mtbl_cf solution_hom_un(source.get_mg(), source.base) ;
105  Mtbl_cf solution_hom_deux(source.get_mg(), source.base) ;
106  Mtbl_cf resultat(source.get_mg(), source.base) ;
107 
108  solution_part.annule_hard() ;
109  solution_hom_un.annule_hard() ;
110  solution_hom_deux.annule_hard() ;
111  resultat.annule_hard() ;
112 
113  // Computation of the SP and SH's in every domain ...
114  int conte = 0 ;
115  for (int zone=0 ; zone<nz ; zone++) {
116 
117  nr = source.get_mg()->get_nr(zone) ;
118  nt = source.get_mg()->get_nt(zone) ;
119  np = source.get_mg()->get_np(zone) ;
120 
121  for (int k=0 ; k<np+1 ; k++)
122  for (int j=0 ; j<nt ; j++) {
123 
124  if (ope_var.operateurs[conte] != 0x0) {
125 
126  // Calcul de la SH
127  sol_hom = new Tbl(ope_var.operateurs[conte]->get_solh()) ;
128 
129  //Calcul de la SP
130  so = new Tbl(nr) ;
131  so->set_etat_qcq() ;
132  for (int i=0 ; i<nr ; i++)
133  so->set(i) = source(zone, k, j, i) ;
134 
135  sol_part = new Tbl(ope_var.operateurs[conte]->get_solp(*so)) ;
136 
137  // Rangement dans les tableaux globaux ;
138  for (int i=0 ; i<nr ; i++) {
139  solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
140  if (sol_hom->get_ndim()==1)
141  solution_hom_un.set(zone, k, j, i) = (*sol_hom)(i) ;
142  else
143  {
144  solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0,i) ;
145  solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1,i) ;
146  }
147  }
148 
149  delete so ;
150  delete sol_hom ;
151  delete sol_part ;
152 
153  }
154  conte ++ ;
155  }
156  }
157 
158 
159  //-------------------------------------------------
160  // ON EST PARTI POUR LE RACCORD (Be carefull ....)
161  //-------------------------------------------------
162 
163  // C'est pas simple toute cette sombre affaire...
164  // Que le cas meme nombre de points dans chaque domaines...
165 
166  int start = 0 ;
167  for (int k=0 ; k<source.get_mg()->get_np(0)+1 ; k++)
168  for (int j=0 ; j<source.get_mg()->get_nt(0) ; j++) {
169  if (ope_var.operateurs[start] != 0x0) {
170 
171  int taille = 2*nz - 2 ;
172  Matrice systeme (taille, taille) ;
173  systeme.set_etat_qcq() ;
174  for (int i=0 ; i<taille ; i++)
175  for (int j2=0 ; j2<taille ; j2++)
176  systeme.set(i,j2) = 0 ;
177  Tbl sec_membre (taille) ;
178  sec_membre.set_etat_qcq() ;
179  for (int i=0 ; i<taille ; i++)
180  sec_membre.set(i) = 0 ;
181 
182  //---------
183  // Noyau :
184  //---------
185  conte = start ;
186 
187  systeme.set(0,0) = ope_var.G_plus(0) *
188  ope_var.operateurs[conte]->val_sh_one_plus() ;
189  systeme.set(1,0) =
190  ope_var.dG_plus(0) * ope_var.operateurs[conte]->val_sh_one_plus() +
191  ope_var.G_plus(0) * ope_var.operateurs[conte]->der_sh_one_plus() ;
192 
193  sec_membre.set(0) -= ope_var.F_plus(0,k,j) +
194  ope_var.G_plus(0) * ope_var.operateurs[conte]->val_sp_plus() ;
195  sec_membre.set(1) -= ope_var.dF_plus(0,k,j) +
196  ope_var.dG_plus(0) * ope_var.operateurs[conte]->val_sp_plus() +
197  ope_var.G_plus(0) * ope_var.operateurs[conte]->der_sp_plus() ;
198 
199  //----------
200  // SHELLS :
201  //----------
202 
203  for (int l=1 ; l<nz-1 ; l++) {
204 
205  // On se met au bon endroit :
206  int np_prec = source.get_mg()->get_np(l-1) ;
207  int nt_prec = source.get_mg()->get_nt(l-1) ;
208  conte += (np_prec+1)*nt_prec ;
209 
210  systeme.set(2*l-2, 2*l-1) = -ope_var.G_minus(l) *
211  ope_var.operateurs[conte]->val_sh_one_minus() ;
212  systeme.set(2*l-2, 2*l) = - ope_var.G_minus(l) *
213  ope_var.operateurs[conte]->val_sh_two_minus() ;
214  systeme.set(2*l-1, 2*l-1) =
215  -ope_var.dG_minus(l)*ope_var.operateurs[conte]->val_sh_one_minus()-
216  ope_var.G_minus(l)*ope_var.operateurs[conte]->der_sh_one_minus() ;
217  systeme.set(2*l-1, 2*l) =
218  -ope_var.dG_minus(l)*ope_var.operateurs[conte]->val_sh_two_minus()-
219  ope_var.G_minus(l)*ope_var.operateurs[conte]->der_sh_two_minus() ;
220 
221  sec_membre.set(2*l-2) += ope_var.F_minus(l,k,j) +
222  ope_var.G_minus(l) * ope_var.operateurs[conte]->val_sp_minus() ;
223  sec_membre.set(2*l-1) += ope_var.dF_minus(l,k,j) +
224  ope_var.dG_minus(l) * ope_var.operateurs[conte]->val_sp_minus() +
225  ope_var.G_minus(l) * ope_var.operateurs[conte]->der_sp_minus() ;
226 
227  // Valeurs en +1 :
228  systeme.set(2*l, 2*l-1) = ope_var.G_plus(l) *
229  ope_var.operateurs[conte]->val_sh_one_plus() ;
230  systeme.set(2*l, 2*l) = ope_var.G_plus(l) *
231  ope_var.operateurs[conte]->val_sh_two_plus() ;
232  systeme.set(2*l+1, 2*l-1) =
233  ope_var.dG_plus(l)*ope_var.operateurs[conte]->val_sh_one_plus()+
234  ope_var.G_plus(l)*ope_var.operateurs[conte]->der_sh_one_plus() ;
235  systeme.set(2*l+1, 2*l) =
236  ope_var.dG_plus(l)*ope_var.operateurs[conte]->val_sh_two_plus()+
237  ope_var.G_plus(l)*ope_var.operateurs[conte]->der_sh_two_plus() ;
238 
239  sec_membre.set(2*l) -= ope_var.F_plus(l,k,j) +
240  ope_var.G_plus(l) * ope_var.operateurs[conte]->val_sp_plus();
241  sec_membre.set(2*l+1) -= ope_var.dF_plus(l,k,j) +
242  ope_var.dG_plus(l) * ope_var.operateurs[conte]->val_sp_plus() +
243  ope_var.G_plus(l) * ope_var.operateurs[conte]->der_sp_plus() ;
244  }
245 
246  //-------
247  // ZEC :
248  //-------
249  int np_prec = source.get_mg()->get_np(nz-2) ;
250  int nt_prec = source.get_mg()->get_nt(nz-2) ;
251  conte += (np_prec+1)*nt_prec ;
252 
253  systeme.set(taille-2, taille-1) = -ope_var.G_minus(nz-1) *
254  ope_var.operateurs[conte]->val_sh_one_minus() ;
255  systeme.set(taille-1, taille-1) =
256  -ope_var.dG_minus(nz-1)*ope_var.operateurs[conte]->val_sh_one_minus()-
257  ope_var.G_minus(nz-1)*ope_var.operateurs[conte]->der_sh_one_minus() ;
258 
259  sec_membre.set(taille-2) += ope_var.F_minus(nz-1,k,j) +
260  ope_var.G_minus(nz-1)*ope_var.operateurs[conte]->val_sp_minus() ;
261  sec_membre.set(taille-1) += ope_var.dF_minus(nz-1,k,j) +
262  ope_var.dG_minus(nz-1) * ope_var.operateurs[conte]->val_sp_minus() +
263  ope_var.G_minus(nz-1) * ope_var.operateurs[conte]->der_sp_minus() ;
264 
265  // On resout le systeme ...
266  if (taille > 2)
267  systeme.set_band(2,2) ;
268  else
269  systeme.set_band(1,1) ;
270 
271  systeme.set_lu() ;
272  Tbl facteur (systeme.inverse(sec_membre)) ;
273 
274  // On range tout ca :
275  // Noyau
276  nr = source.get_mg()->get_nr(0) ;
277  for (int i=0 ; i<nr ; i++)
278  resultat.set(0,k,j,i) = solution_part(0,k,j,i) +facteur(0)*solution_hom_un(0,k,j,i) ;
279 
280  // Shells
281  for (int l=1 ; l<nz-1 ; l++) {
282  nr = source.get_mg()->get_nr(l) ;
283  for (int i=0 ; i<nr ; i++)
284  resultat.set(l,k,j,i) = solution_part(l,k,j,i) +
285  facteur(2*l-1)*solution_hom_un(l,k,j,i) +
286  facteur(2*l)*solution_hom_deux(l,k,j,i) ;
287  }
288 
289  // Zec
290  nr = source.get_mg()->get_nr(nz-1) ;
291  for (int i=0 ; i<nr ; i++)
292  resultat.set(nz-1,k,j,i) = solution_part(nz-1,k,j,i) +
293  facteur(taille-1)*solution_hom_un(nz-1,k,j,i) ;
294  }
295  start ++ ;
296  }
297 
298  return resultat;
299 }
300 
301 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Mtbl_cf::get_mg
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition: mtbl_cf.h:453