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 
21 char 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 
72 namespace Lorene {
73 Mtbl_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 }
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