LORENE
sol_elliptic_fixe_der_zero.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_fixe_der_zeroC[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/sol_elliptic_fixe_der_zero.C,v 1.5 2014/10/13 08:53:30 j_novak Exp $" ;
22 
23 /*
24  * $Id: sol_elliptic_fixe_der_zero.C,v 1.5 2014/10/13 08:53:30 j_novak Exp $
25  *
26  */
27 
28 // Header C :
29 #include <cstdlib>
30 #include <cmath>
31 
32 // Headers Lorene :
33 #include "tbl.h"
34 #include "mtbl_cf.h"
35 #include "map.h"
36 #include "param_elliptic.h"
37 
38 
39  //----------------------------------------------
40  // Version Mtbl_cf
41  //----------------------------------------------
42 
43 
44 
45 namespace Lorene {
46 Mtbl_cf elliptic_solver_fixe_der_zero (double valeur,
47  const Param_elliptic& ope_var,
48  const Mtbl_cf& source) {
49  // Verifications d'usage sur les zones
50  int nz = source.get_mg()->get_nzone() ;
51  assert (nz>1) ;
52  assert (source.get_mg()->get_type_r(0) == RARE) ;
53  assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
54  for (int l=1 ; l<nz-1 ; l++)
55  assert(source.get_mg()->get_type_r(l) == FIN) ;
56 
57  // donnees sur la zone
58  int nr, nt, np ;
59 
60  //Rangement des valeurs intermediaires
61  Tbl *so ;
62  Tbl *sol_hom ;
63  Tbl *sol_part ;
64 
65 
66  // Rangement des solutions, avant raccordement
67  Mtbl_cf solution_part(source.get_mg(), source.base) ;
68  Mtbl_cf solution_hom_un(source.get_mg(), source.base) ;
69  Mtbl_cf solution_hom_deux(source.get_mg(), source.base) ;
70  Mtbl_cf resultat(source.get_mg(), source.base) ;
71 
72  solution_part.annule_hard() ;
73  solution_hom_un.annule_hard() ;
74  solution_hom_deux.annule_hard() ;
75  resultat.annule_hard() ;
76 
77  // Computation of the SP and SH's in every domain ...
78  int conte = 0 ;
79  for (int zone=0 ; zone<nz ; zone++) {
80  nr = source.get_mg()->get_nr(zone) ;
81  nt = source.get_mg()->get_nt(zone) ;
82  np = source.get_mg()->get_np(zone) ;
83 
84  for (int k=0 ; k<np+1 ; k++)
85  for (int j=0 ; j<nt ; j++) {
86  if (ope_var.operateurs[conte] != 0x0) {
87 
88  // Calcul de la SH
89  sol_hom = new Tbl(ope_var.operateurs[conte]->get_solh()) ;
90 
91  //Calcul de la SP
92  so = new Tbl(nr) ;
93  so->set_etat_qcq() ;
94  for (int i=0 ; i<nr ; i++)
95  so->set(i) = source(zone, k, j, i) ;
96 
97  sol_part = new Tbl(ope_var.operateurs[conte]->get_solp(*so)) ;
98 
99  // Rangement dans les tableaux globaux ;
100  for (int i=0 ; i<nr ; i++) {
101  solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
102  if (sol_hom->get_ndim()==1)
103  solution_hom_un.set(zone, k, j, i) = (*sol_hom)(i) ;
104  else
105  {
106  solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0,i) ;
107  solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1,i) ;
108  }
109  }
110 
111  delete so ;
112  delete sol_hom ;
113  delete sol_part ;
114 
115  }
116  conte ++ ;
117  }
118  }
119 
120  //-------------------------------------------------
121  // ON EST PARTI POUR LE RACCORD (Be carefull ....)
122  //-------------------------------------------------
123 
124  // C'est pas simple toute cette sombre affaire...
125  // POUR LE MOMENT QUE LE CAS l==0 ;
126  // Que le cas meme nombre de points dans chaque domaines...
127 
128  int start = 0 ;
129  for (int k=0 ; k<source.get_mg()->get_np(0)+1 ; k++)
130  for (int j=0 ; j<source.get_mg()->get_nt(0) ; j++) {
131  if (ope_var.operateurs[start] != 0x0) {
132 
133  int taille = 2*nz - 2 ;
134  Matrice systeme (taille, taille) ;
135  systeme.set_etat_qcq() ;
136  for (int i=0 ; i<taille ; i++)
137  for (int j2=0 ; j2<taille ; j2++)
138  systeme.set(i,j2) = 0 ;
139  Tbl sec_membre (taille) ;
140  sec_membre.set_etat_qcq() ;
141  for (int i=0 ; i<taille ; i++)
142  sec_membre.set(i) = 0 ;
143 
144  //---------
145  // Noyau :
146  //---------
147  conte = start ;
148 
149  systeme.set(0,0) = ope_var.G_plus(0) *
150  ope_var.operateurs[conte]->val_sh_one_plus() ;
151 
152  // On relache derivee
153  systeme.set(1,0) =
154  ope_var.dG_minus(0) * ope_var.operateurs[conte]->val_sh_one_minus() +
155  ope_var.G_minus(0) * ope_var.operateurs[conte]->der_sh_one_minus() ;
156 
157  sec_membre.set(0) -= ope_var.F_plus(0,k,j) +
158  ope_var.G_plus(0) * ope_var.operateurs[conte]->val_sp_plus() ;
159 
160  if ((k==0) && (j==0))
161  sec_membre.set(1) -= -valeur +
162  ope_var.dF_minus(0,k,j) +
163  ope_var.dG_minus(0) * ope_var.operateurs[conte]->val_sp_minus() +
164  ope_var.G_minus(0) * ope_var.operateurs[conte]->der_sp_minus() ;
165 
166  //----------
167  // SHELLS :
168  //----------
169 
170  for (int l=1 ; l<nz-1 ; l++) {
171 
172  // On se met au bon endroit :
173  int np_prec = source.get_mg()->get_np(l-1) ;
174  int nt_prec = source.get_mg()->get_nt(l-1) ;
175  conte += (np_prec+1)*nt_prec ;
176 
177 
178  systeme.set(2*l-2, 2*l-1) = -ope_var.G_minus(l) *
179  ope_var.operateurs[conte]->val_sh_one_minus() ;
180  systeme.set(2*l-2, 2*l) = - ope_var.G_minus(l) *
181  ope_var.operateurs[conte]->val_sh_two_minus() ;
182  if ((l!=1) || (k!=0) || (j!=0)) {
183  systeme.set(2*l-1, 2*l-1) =
184  -ope_var.dG_minus(l)*ope_var.operateurs[conte]->val_sh_one_minus()-
185  ope_var.G_minus(l)*ope_var.operateurs[conte]->der_sh_one_minus() ;
186  systeme.set(2*l-1, 2*l) =
187  -ope_var.dG_minus(l)*ope_var.operateurs[conte]->val_sh_two_minus()-
188  ope_var.G_minus(l)*ope_var.operateurs[conte]->der_sh_two_minus() ;
189  }
190  sec_membre.set(2*l-2) += ope_var.F_minus(l,k,j) +
191  ope_var.G_minus(l) * ope_var.operateurs[conte]->val_sp_minus() ;
192  if ((l!=1) || (k!=0) || (j!=0)) {
193  sec_membre.set(2*l-1) += ope_var.dF_minus(l,k,j) +
194  ope_var.dG_minus(l) * ope_var.operateurs[conte]->val_sp_minus() +
195  ope_var.G_minus(l) * ope_var.operateurs[conte]->der_sp_minus() ;
196  }
197 
198  // Valeurs en +1 :
199 
200  systeme.set(2*l, 2*l-1) = ope_var.G_plus(l) *
201  ope_var.operateurs[conte]->val_sh_one_plus() ;
202  systeme.set(2*l, 2*l) = ope_var.G_plus(l) *
203  ope_var.operateurs[conte]->val_sh_two_plus() ;
204  systeme.set(2*l+1, 2*l-1) =
205  ope_var.dG_plus(l)*ope_var.operateurs[conte]->val_sh_one_plus()+
206  ope_var.G_plus(l)*ope_var.operateurs[conte]->der_sh_one_plus() ;
207  systeme.set(2*l+1, 2*l) =
208  ope_var.dG_plus(l)*ope_var.operateurs[conte]->val_sh_two_plus()+
209  ope_var.G_plus(l)*ope_var.operateurs[conte]->der_sh_two_plus() ;
210 
211  sec_membre.set(2*l) -= ope_var.F_plus(l,k,j) +
212  ope_var.G_plus(l) * ope_var.operateurs[conte]->val_sp_plus();
213  sec_membre.set(2*l+1) -= ope_var.dF_plus(l,k,j) +
214  ope_var.dG_plus(l) * ope_var.operateurs[conte]->val_sp_plus() +
215  ope_var.G_plus(l) * ope_var.operateurs[conte]->der_sp_plus() ;
216  }
217 
218  //-------
219  // ZEC :
220  //-------
221  int np_prec = source.get_mg()->get_np(nz-2) ;
222  int nt_prec = source.get_mg()->get_nt(nz-2) ;
223  conte += (np_prec+1)*nt_prec ;
224 
225  systeme.set(taille-2, taille-1) = -ope_var.G_minus(nz-1) *
226  ope_var.operateurs[conte]->val_sh_one_minus() ;
227  systeme.set(taille-1, taille-1) =
228  -ope_var.dG_minus(nz-1)*ope_var.operateurs[conte]->val_sh_one_minus()-
229  ope_var.G_minus(nz-1)*ope_var.operateurs[conte]->der_sh_one_minus() ;
230 
231  sec_membre.set(taille-2) += ope_var.F_minus(nz-1,k,j) +
232  ope_var.G_minus(nz-1)*ope_var.operateurs[conte]->val_sp_minus() ;
233  sec_membre.set(taille-1) += ope_var.dF_minus(nz-1,k,j) +
234  ope_var.dG_minus(nz-1) * ope_var.operateurs[conte]->val_sp_minus() +
235  ope_var.G_minus(nz-1) * ope_var.operateurs[conte]->der_sp_minus() ;
236 
237  // On resout le systeme ...
238  if (taille > 2)
239  systeme.set_band(2,2) ;
240  else
241  systeme.set_band(1,1) ;
242 
243  systeme.set_lu() ;
244  Tbl facteur (systeme.inverse(sec_membre)) ;
245 
246  // On range tout ca :
247  // Noyau
248  nr = source.get_mg()->get_nr(0) ;
249  for (int i=0 ; i<nr ; i++)
250  resultat.set(0,k,j,i) = solution_part(0,k,j,i)
251  +facteur(0)*solution_hom_un(0,k,j,i) ;
252 
253  // Shells
254  for (int l=1 ; l<nz-1 ; l++) {
255  nr = source.get_mg()->get_nr(l) ;
256  for (int i=0 ; i<nr ; i++)
257  resultat.set(l,k,j,i) = solution_part(l,k,j,i) +
258  facteur(2*l-1)*solution_hom_un(l,k,j,i) +
259  facteur(2*l)*solution_hom_deux(l,k,j,i) ;
260  }
261 
262  // Zec
263  nr = source.get_mg()->get_nr(nz-1) ;
264  for (int i=0 ; i<nr ; i++)
265  resultat.set(nz-1,k,j,i) = solution_part(nz-1,k,j,i) +
266  facteur(taille-1)*solution_hom_un(nz-1,k,j,i) ;
267  }
268  start ++ ;
269  }
270 
271  return resultat;
272 }
273 
274 }
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition: mtbl_cf.h:463
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465