LORENE
sol_elliptic_boundary.C
1 /*
2  * Copyright (c) 2003 Philippe Grandclement
3  * Jose Luis Jaramillo
4  * Francois Limousin
5  *
6  * This file is part of LORENE.
7  *
8  * LORENE is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License version 2
10  * as published by the Free Software Foundation.
11  *
12  * LORENE is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with LORENE; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20  *
21  */
22 
23 char sol_elliptic_boundary_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/sol_elliptic_boundary.C,v 1.4 2014/10/13 08:53:30 j_novak Exp $" ;
24 
25 /*
26  * $Id: sol_elliptic_boundary.C,v 1.4 2014/10/13 08:53:30 j_novak Exp $
27  * $Log: sol_elliptic_boundary.C,v $
28  * Revision 1.4 2014/10/13 08:53:30 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.3 2014/10/06 15:16:10 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.2 2008/08/20 15:03:55 n_vasset
35  * Correction on how the boundary condition is imposed
36  *
37  * Revision 1.1 2005/06/09 08:01:05 f_limousin
38  * Implement a new function sol_elliptic_boundary() and
39  * Vector::poisson_boundary(...) which solve the vectorial poisson
40  * equation (method 6) with an inner boundary condition.
41  *
42  *
43  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/sol_elliptic_boundary.C,v 1.4 2014/10/13 08:53:30 j_novak Exp $
44  *
45  */
46 
47 // Header C :
48 #include <cstdlib>
49 #include <cmath>
50 
51 // Headers Lorene :
52 #include "param_elliptic.h"
53 #include "tbl.h"
54 #include "mtbl_cf.h"
55 #include "map.h"
56 
57 
58  //----------------------------------------------
59  // Version Mtbl_cf
60  //----------------------------------------------
61 
62 
63 
64 namespace Lorene {
65 Mtbl_cf elliptic_solver_boundary (const Param_elliptic& ope_var, const Mtbl_cf& source,
66  const Mtbl_cf& bound, double fact_dir, double fact_neu ) {
67  // Verifications d'usage sur les zones
68  int nz = source.get_mg()->get_nzone() ;
69  assert (nz>1) ;
70  assert (source.get_mg()->get_type_r(0) == RARE) ;
71  assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
72  for (int l=1 ; l<nz-1 ; l++)
73  assert(source.get_mg()->get_type_r(l) == FIN) ;
74 
75  // donnees sur la zone
76  int nr, nt, np ;
77 
78  //Rangement des valeurs intermediaires
79  Tbl *so ;
80  Tbl *sol_hom ;
81  Tbl *sol_part ;
82 
83 
84  // Rangement des solutions, avant raccordement
85  Mtbl_cf solution_part(source.get_mg(), source.base) ;
86  Mtbl_cf solution_hom_un(source.get_mg(), source.base) ;
87  Mtbl_cf solution_hom_deux(source.get_mg(), source.base) ;
88  Mtbl_cf resultat(source.get_mg(), source.base) ;
89 
90  solution_part.annule_hard() ;
91  solution_hom_un.annule_hard() ;
92  solution_hom_deux.annule_hard() ;
93  resultat.annule_hard() ;
94 
95  // Computation of the SP and SH's in every domain ...
96  int conte = 0 ;
97  for (int zone=0 ; zone<nz ; zone++) {
98 
99  nr = source.get_mg()->get_nr(zone) ;
100  nt = source.get_mg()->get_nt(zone) ;
101  np = source.get_mg()->get_np(zone) ;
102 
103  for (int k=0 ; k<np+1 ; k++)
104  for (int j=0 ; j<nt ; j++) {
105 
106  if (ope_var.operateurs[conte] != 0x0) {
107 
108  // Calcul de la SH
109  sol_hom = new Tbl(ope_var.operateurs[conte]->get_solh()) ;
110 
111  //Calcul de la SP
112  so = new Tbl(nr) ;
113  so->set_etat_qcq() ;
114  for (int i=0 ; i<nr ; i++)
115  so->set(i) = source(zone, k, j, i) ;
116 
117  sol_part = new Tbl(ope_var.operateurs[conte]->get_solp(*so)) ;
118 
119  // Rangement dans les tableaux globaux ;
120  for (int i=0 ; i<nr ; i++) {
121  solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
122  if (sol_hom->get_ndim()==1)
123  solution_hom_un.set(zone, k, j, i) = (*sol_hom)(i) ;
124  else
125  {
126  solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0,i) ;
127  solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1,i) ;
128  }
129  }
130 
131  delete so ;
132  delete sol_hom ;
133  delete sol_part ;
134 
135  }
136  conte ++ ;
137  }
138  }
139 
140 
141  //-------------------------------------------------
142  // ON EST PARTI POUR LE RACCORD (Be careful ....)
143  //-------------------------------------------------
144 
145  // C'est pas simple toute cette sombre affaire...
146  // Que de cas meme nombre de points dans chaque domaines...
147 
148  int start = 0 ;
149  for (int k=0 ; k<source.get_mg()->get_np(0)+1 ; k++)
150  for (int j=0 ; j<source.get_mg()->get_nt(0) ; j++) {
151  if (ope_var.operateurs[start] != 0x0) {
152 
153  int taille = 2*nz - 3 ;
154  Matrice systeme (taille, taille) ;
155  systeme.set_etat_qcq() ;
156  for (int i=0 ; i<taille ; i++)
157  for (int j2=0 ; j2<taille ; j2++)
158  systeme.set(i,j2) = 0 ;
159  Tbl sec_membre (taille) ;
160  sec_membre.set_etat_qcq() ;
161  for (int i=0 ; i<taille ; i++)
162  sec_membre.set(i) = 0 ;
163 
164  //----------
165  // BOUNDARY :
166  //----------
167  conte = start ;
168 
169  // Boundary value at an angular point
170 
171  // Setting the right hand side term
172 
173  sec_membre.set(0) -= bound.val_in_bound_jk(1, j, k)/sqrt(2.) ; //Relevant value is in the first shell and only in the first coefficient
174 
175 
176  //--------------
177  // FIRST SHELL :
178  //--------------
179 
180  int l_1=1 ;
181 
182  // On se met au bon endroit :
183  int np_prec_1 = source.get_mg()->get_np(l_1-1) ;
184  int nt_prec_1 = source.get_mg()->get_nt(l_1-1) ;
185  conte += (np_prec_1+1)*nt_prec_1 ;
186 
187  systeme.set(0, 0) = fact_dir * (-ope_var.G_minus(l_1) *
188  ope_var.operateurs[conte]->val_sh_one_minus() )
189  + fact_neu *
190  ( -ope_var.dG_minus(l_1)*ope_var.operateurs[conte]->val_sh_one_minus()-
191  ope_var.G_minus(l_1)*ope_var.operateurs[conte]->der_sh_one_minus() );
192 
193  systeme.set(0, 1) = fact_dir * (- ope_var.G_minus(l_1) *
194  ope_var.operateurs[conte]->val_sh_two_minus() ) + fact_neu *
195  (-ope_var.dG_minus(l_1)*ope_var.operateurs[conte]->val_sh_two_minus()-
196  ope_var.G_minus(l_1)*ope_var.operateurs[conte]->der_sh_two_minus() ) ;
197 
198  //Completing the right hand side term
199  sec_membre.set(0) += fact_dir * (ope_var.F_minus(l_1,k,j) +
200  ope_var.G_minus(l_1) * ope_var.operateurs[conte]->val_sp_minus() ) +
201  fact_neu * ( ope_var.dF_minus(l_1,k,j) +
202  ope_var.dG_minus(l_1) * ope_var.operateurs[conte]->val_sp_minus() +
203  ope_var.G_minus(l_1) * ope_var.operateurs[conte]->der_sp_minus() ) ;
204 
205 
206  // Valeurs en +1 :
207  systeme.set(2*l_1-1, 2*l_1-2) = ope_var.G_plus(l_1) *
208  ope_var.operateurs[conte]->val_sh_one_plus() ;
209  systeme.set(2*l_1-1, 2*l_1-1) = ope_var.G_plus(l_1) *
210  ope_var.operateurs[conte]->val_sh_two_plus() ;
211  systeme.set(2*l_1, 2*l_1-2) =
212  ope_var.dG_plus(l_1)*ope_var.operateurs[conte]->val_sh_one_plus()+
213  ope_var.G_plus(l_1)*ope_var.operateurs[conte]->der_sh_one_plus() ;
214  systeme.set(2*l_1, 2*l_1-1) =
215  ope_var.dG_plus(l_1)*ope_var.operateurs[conte]->val_sh_two_plus()+
216  ope_var.G_plus(l_1)*ope_var.operateurs[conte]->der_sh_two_plus() ;
217 
218  sec_membre.set(2*l_1-1) -= ope_var.F_plus(l_1,k,j) +
219  ope_var.G_plus(l_1) * ope_var.operateurs[conte]->val_sp_plus();
220  sec_membre.set(2*l_1) -= ope_var.dF_plus(l_1,k,j) +
221  ope_var.dG_plus(l_1) * ope_var.operateurs[conte]->val_sp_plus() +
222  ope_var.G_plus(l_1) * ope_var.operateurs[conte]->der_sp_plus() ;
223 
224 
225 
226  //----------
227  // SHELLS :
228  //----------
229 // assert (nz-1 > 2) ; // At least two shells
230  for (int l=2 ; l<nz-1 ; l++) {
231 
232  // On se met au bon endroit :
233  int np_prec = source.get_mg()->get_np(l-1) ;
234  int nt_prec = source.get_mg()->get_nt(l-1) ;
235  conte += (np_prec+1)*nt_prec ;
236 
237  systeme.set(2*l-3, 2*l-2) = -ope_var.G_minus(l) *
238  ope_var.operateurs[conte]->val_sh_one_minus() ;
239  systeme.set(2*l-3, 2*l-1) = - ope_var.G_minus(l) *
240  ope_var.operateurs[conte]->val_sh_two_minus() ;
241  systeme.set(2*l-2, 2*l-2) =
242  -ope_var.dG_minus(l)*ope_var.operateurs[conte]->val_sh_one_minus()-
243  ope_var.G_minus(l)*ope_var.operateurs[conte]->der_sh_one_minus() ;
244  systeme.set(2*l-2, 2*l-1) =
245  -ope_var.dG_minus(l)*ope_var.operateurs[conte]->val_sh_two_minus()-
246  ope_var.G_minus(l)*ope_var.operateurs[conte]->der_sh_two_minus() ;
247 
248  sec_membre.set(2*l-3) += ope_var.F_minus(l,k,j) +
249  ope_var.G_minus(l) * ope_var.operateurs[conte]->val_sp_minus() ;
250  sec_membre.set(2*l-2) += ope_var.dF_minus(l,k,j) +
251  ope_var.dG_minus(l) * ope_var.operateurs[conte]->val_sp_minus() +
252  ope_var.G_minus(l) * ope_var.operateurs[conte]->der_sp_minus() ;
253 
254  // Valeurs en +1 :
255  systeme.set(2*l-1, 2*l-2) = ope_var.G_plus(l) *
256  ope_var.operateurs[conte]->val_sh_one_plus() ;
257  systeme.set(2*l-1, 2*l-1) = ope_var.G_plus(l) *
258  ope_var.operateurs[conte]->val_sh_two_plus() ;
259  systeme.set(2*l, 2*l-2) =
260  ope_var.dG_plus(l)*ope_var.operateurs[conte]->val_sh_one_plus()+
261  ope_var.G_plus(l)*ope_var.operateurs[conte]->der_sh_one_plus() ;
262  systeme.set(2*l, 2*l-1) =
263  ope_var.dG_plus(l)*ope_var.operateurs[conte]->val_sh_two_plus()+
264  ope_var.G_plus(l)*ope_var.operateurs[conte]->der_sh_two_plus() ;
265 
266  sec_membre.set(2*l-1) -= ope_var.F_plus(l,k,j) +
267  ope_var.G_plus(l) * ope_var.operateurs[conte]->val_sp_plus();
268  sec_membre.set(2*l) -= ope_var.dF_plus(l,k,j) +
269  ope_var.dG_plus(l) * ope_var.operateurs[conte]->val_sp_plus() +
270  ope_var.G_plus(l) * ope_var.operateurs[conte]->der_sp_plus() ;
271  }
272 
273  //-------
274  // ZEC :
275  //-------
276  int np_prec = source.get_mg()->get_np(nz-2) ;
277  int nt_prec = source.get_mg()->get_nt(nz-2) ;
278  conte += (np_prec+1)*nt_prec ;
279 
280  systeme.set(taille-2, taille-1) = -ope_var.G_minus(nz-1) *
281  ope_var.operateurs[conte]->val_sh_one_minus() ;
282  systeme.set(taille-1, taille-1) =
283  -ope_var.dG_minus(nz-1)*ope_var.operateurs[conte]->val_sh_one_minus()-
284  ope_var.G_minus(nz-1)*ope_var.operateurs[conte]->der_sh_one_minus() ;
285 
286  sec_membre.set(taille-2) += ope_var.F_minus(nz-1,k,j) +
287  ope_var.G_minus(nz-1)*ope_var.operateurs[conte]->val_sp_minus() ;
288  sec_membre.set(taille-1) += ope_var.dF_minus(nz-1,k,j) +
289  ope_var.dG_minus(nz-1) * ope_var.operateurs[conte]->val_sp_minus() +
290  ope_var.G_minus(nz-1) * ope_var.operateurs[conte]->der_sp_minus() ;
291 
292  // On resout le systeme ...
293  if (taille > 2)
294  systeme.set_band(2,2) ;
295  else
296  systeme.set_band(1,1) ;
297 
298  systeme.set_lu() ;
299  Tbl facteur (systeme.inverse(sec_membre)) ;
300 
301  // On range tout ca :
302 
303  // Shells
304  for (int l=1 ; l<nz-1 ; l++) {
305  nr = source.get_mg()->get_nr(l) ;
306  for (int i=0 ; i<nr ; i++)
307  resultat.set(l,k,j,i) = solution_part(l,k,j,i) +
308  facteur(2*l-2)*solution_hom_un(l,k,j,i) +
309  facteur(2*l-1)*solution_hom_deux(l,k,j,i) ;
310  }
311 
312  // Zec
313  nr = source.get_mg()->get_nr(nz-1) ;
314  for (int i=0 ; i<nr ; i++)
315  resultat.set(nz-1,k,j,i) = solution_part(nz-1,k,j,i) +
316  facteur(taille-1)*solution_hom_un(nz-1,k,j,i) ;
317  }
318  start ++ ;
319  }
320 
321  return resultat;
322 }
323 
324 }
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene prototypes.
Definition: app_hor.h:64
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition: mtbl_cf.h:453
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448