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