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
21char 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
45namespace Lorene {
46Mtbl_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}
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
Lorene prototypes.
Definition app_hor.h:64