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
21char 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
83namespace Lorene {
84Mtbl_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}
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition mtbl_cf.h:453
Lorene prototypes.
Definition app_hor.h:64