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
23char 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
64namespace Lorene {
65Mtbl_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}
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
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene prototypes.
Definition app_hor.h:64