57char sol_elliptic_sin_zec_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/sol_elliptic_sin_zec.C,v 1.8 2014/10/13 08:53:30 j_novak Exp $" ;
62#include <gsl/gsl_sf_bessel.h>
68#include "param_elliptic.h"
76Mtbl_cf elliptic_solver_sin_zec (
const Param_elliptic& ope_var,
77 const Mtbl_cf& source,
double* amplis,
double* phases) {
82 assert (source.get_mg()->get_type_r(0) == RARE) ;
83 assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
84 for (
int l=1 ; l<nz-1 ; l++)
85 assert(source.get_mg()->get_type_r(l) == FIN) ;
97 Mtbl_cf solution_part(source.get_mg(), source.base) ;
98 Mtbl_cf solution_hom_un(source.get_mg(), source.base) ;
99 Mtbl_cf solution_hom_deux(source.get_mg(), source.base) ;
100 Mtbl_cf resultat(source.get_mg(), source.base) ;
102 solution_part.annule_hard() ;
103 solution_hom_un.annule_hard() ;
104 solution_hom_deux.annule_hard() ;
105 resultat.annule_hard() ;
109 for (
int zone=0 ; zone<nz-1 ; zone++) {
110 nr = source.get_mg()->get_nr(zone) ;
111 nt = source.get_mg()->get_nt(zone) ;
112 np = source.get_mg()->get_np(zone) ;
114 for (
int k=0 ; k<np+1 ; k++)
115 for (
int j=0 ; j<nt ; j++) {
116 if (ope_var.operateurs[conte] != 0x0) {
118 sol_hom =
new Tbl(ope_var.operateurs[conte]->get_solh()) ;
123 for (
int i=0 ; i<nr ; i++)
124 so->set(i) = source(zone, k, j, i) ;
126 sol_part =
new Tbl(ope_var.operateurs[conte]->get_solp(*so)) ;
129 for (
int i=0 ; i<nr ; i++) {
130 solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
131 if (sol_hom->get_ndim()==1)
132 solution_hom_un.set(zone, k, j, i) = (*sol_hom)(i) ;
135 solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0,i) ;
136 solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1,i) ;
157 for (
int k=0 ; k< source.get_mg()->get_np(0)+1; k++)
158 for (
int j=0 ; j<source.get_mg()->get_nt(0) ; j++) {
159 if (ope_var.operateurs[start] != 0x0) {
161 int taille = 2*nz - 2 ;
162 Matrice systeme (taille, taille) ;
163 systeme.set_etat_qcq() ;
164 for (
int i=0 ; i<taille ; i++)
165 for (
int j2=0 ; j2<taille ; j2++)
166 systeme.set(i,j2) = 0 ;
167 Tbl sec_membre (taille) ;
168 sec_membre.set_etat_qcq() ;
169 for (
int i=0 ; i<taille ; i++)
170 sec_membre.set(i) = 0 ;
176 systeme.set(0,0) = ope_var.G_plus(0) *
177 ope_var.operateurs[conte]->val_sh_one_plus() ;
179 ope_var.dG_plus(0) * ope_var.operateurs[conte]->val_sh_one_plus() +
180 ope_var.G_plus(0) * ope_var.operateurs[conte]->der_sh_one_plus() ;
182 sec_membre.set(0) -= ope_var.F_plus(0,k,j) +
183 ope_var.G_plus(0) * ope_var.operateurs[conte]->val_sp_plus() ;
184 sec_membre.set(1) -= ope_var.dF_plus(0,k,j) +
185 ope_var.dG_plus(0) * ope_var.operateurs[conte]->val_sp_plus() +
186 ope_var.G_plus(0) * ope_var.operateurs[conte]->der_sp_plus() ;
192 for (
int l=1 ; l<nz-1 ; l++) {
195 int np_prec = source.get_mg()->get_np(l-1) ;
196 int nt_prec = source.get_mg()->get_nt(l-1) ;
197 conte += (np_prec+1)*nt_prec ;
199 systeme.set(2*l-2, 2*l-1) = -ope_var.G_minus(l) *
200 ope_var.operateurs[conte]->val_sh_one_minus() ;
201 systeme.set(2*l-2, 2*l) = - ope_var.G_minus(l) *
202 ope_var.operateurs[conte]->val_sh_two_minus() ;
203 systeme.set(2*l-1, 2*l-1) =
204 -ope_var.dG_minus(l)*ope_var.operateurs[conte]->val_sh_one_minus()-
205 ope_var.G_minus(l)*ope_var.operateurs[conte]->der_sh_one_minus() ;
206 systeme.set(2*l-1, 2*l) =
207 -ope_var.dG_minus(l)*ope_var.operateurs[conte]->val_sh_two_minus()-
208 ope_var.G_minus(l)*ope_var.operateurs[conte]->der_sh_two_minus() ;
210 sec_membre.set(2*l-2) += ope_var.F_minus(l,k,j) +
211 ope_var.G_minus(l) * ope_var.operateurs[conte]->val_sp_minus() ;
212 sec_membre.set(2*l-1) += ope_var.dF_minus(l,k,j) +
213 ope_var.dG_minus(l) * ope_var.operateurs[conte]->val_sp_minus() +
214 ope_var.G_minus(l) * ope_var.operateurs[conte]->der_sp_minus() ;
217 systeme.set(2*l, 2*l-1) = ope_var.G_plus(l) *
218 ope_var.operateurs[conte]->val_sh_one_plus() ;
219 systeme.set(2*l, 2*l) = ope_var.G_plus(l) *
220 ope_var.operateurs[conte]->val_sh_two_plus() ;
222 systeme.set(2*l+1, 2*l-1) =
223 ope_var.dG_plus(l)*ope_var.operateurs[conte]->val_sh_one_plus()+
224 ope_var.G_plus(l)*ope_var.operateurs[conte]->der_sh_one_plus() ;
225 systeme.set(2*l+1, 2*l) =
226 ope_var.dG_plus(l)*ope_var.operateurs[conte]->val_sh_two_plus()+
227 ope_var.G_plus(l)*ope_var.operateurs[conte]->der_sh_two_plus() ;
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();
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() ;
239 double val_sh =
cos(phases[start])*ope_var.operateurs[conte]->val_sh_one_plus()
240 +
sin(phases[start])*ope_var.operateurs[conte]->val_sh_two_plus() ;
241 double der_sh =
cos(phases[start])*ope_var.operateurs[conte]->der_sh_one_plus()
242 +
sin(phases[start])*ope_var.operateurs[conte]->der_sh_two_plus() ;
244 systeme.
set(taille-2, taille-1) = -ope_var.G_minus(nz-1) * val_sh ;
245 systeme.
set(taille-1, taille-1) = -ope_var.dG_minus(nz-1)*val_sh- ope_var.G_minus(nz-1)*der_sh ;
247 sec_membre.set(taille-2) += ope_var.F_minus(nz-1,k,j) ;
248 sec_membre.set(taille-1) += ope_var.dF_minus(nz-1,k,j) ;
252 systeme.set_band(2,2) ;
254 systeme.set_band(1,1) ;
257 Tbl facteur (systeme.inverse(sec_membre)) ;
259 amplis[start] = facteur(taille-1) ;
263 nr = source.get_mg()->get_nr(0) ;
264 for (
int i=0 ; i<nr ; i++)
265 resultat.set(0,k,j,i) = solution_part(0,k,j,i)
266 +facteur(0)*solution_hom_un(0,k,j,i) ;
269 for (
int l=1 ; l<nz-1 ; l++) {
270 nr = source.get_mg()->get_nr(l) ;
271 for (
int i=0 ; i<nr ; i++)
272 resultat.set(l,k,j,i) = solution_part(l,k,j,i) +
273 facteur(2*l-1)*solution_hom_un(l,k,j,i) +
274 facteur(2*l)*solution_hom_deux(l,k,j,i) ;
Tbl & set(int l)
Read/write of the value in a given domain.
int get_nzone() const
Returns the number of domains.
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Cmp sin(const Cmp &)
Sine.
Cmp cos(const Cmp &)
Cosine.