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 $" ;
52#include "param_elliptic.h"
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 ) {
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) ;
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) ;
90 solution_part.annule_hard() ;
91 solution_hom_un.annule_hard() ;
92 solution_hom_deux.annule_hard() ;
93 resultat.annule_hard() ;
97 for (
int zone=0 ; zone<nz ; zone++) {
99 nr = source.get_mg()->get_nr(zone) ;
100 nt = source.get_mg()->get_nt(zone) ;
101 np = source.get_mg()->get_np(zone) ;
103 for (
int k=0 ; k<np+1 ; k++)
104 for (
int j=0 ; j<nt ; j++) {
106 if (ope_var.operateurs[conte] != 0x0) {
109 sol_hom =
new Tbl(ope_var.operateurs[conte]->get_solh()) ;
114 for (
int i=0 ; i<nr ; i++)
115 so->set(i) = source(zone, k, j, i) ;
117 sol_part =
new Tbl(ope_var.operateurs[conte]->get_solp(*so)) ;
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) ;
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) ;
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) {
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 ;
173 sec_membre.set(0) -= bound.val_in_bound_jk(1, j, k)/
sqrt(2.) ;
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 ;
187 systeme.set(0, 0) = fact_dir * (-ope_var.G_minus(l_1) *
188 ope_var.operateurs[conte]->val_sh_one_minus() )
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() );
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() ) ;
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() ) ;
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() ;
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() ;
230 for (
int l=2 ; l<nz-1 ; l++) {
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 ;
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() ;
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() ;
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() ;
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() ;
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 ;
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() ;
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() ;
294 systeme.set_band(2,2) ;
296 systeme.set_band(1,1) ;
299 Tbl facteur (systeme.inverse(sec_membre)) ;
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) ;
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) ;
int get_nzone() const
Returns the number of domains.
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Cmp sqrt(const Cmp &)
Square root.