LORENE
poisson_frontiere_double.C
1/*
2 * Copyright (c) 2000-2001 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 as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23char poisson_frontiere_double_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_frontiere_double.C,v 1.3 2014/10/13 08:53:29 j_novak Exp $" ;
24
25/*
26 * $Id: poisson_frontiere_double.C,v 1.3 2014/10/13 08:53:29 j_novak Exp $
27 * $Log: poisson_frontiere_double.C,v $
28 * Revision 1.3 2014/10/13 08:53:29 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.2 2014/10/06 15:16:09 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
35 * LORENE
36 *
37 * Revision 2.1 2000/05/15 15:46:43 phil
38 * *** empty log message ***
39 *
40 * Revision 2.0 2000/04/27 15:19:52 phil
41 * *** empty log message ***
42 *
43 *
44 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_frontiere_double.C,v 1.3 2014/10/13 08:53:29 j_novak Exp $
45 *
46 */
47
48
49// Header C :
50#include <cstdlib>
51#include <cmath>
52
53// Headers Lorene :
54#include "matrice.h"
55#include "tbl.h"
56#include "mtbl_cf.h"
57#include "map.h"
58#include "base_val.h"
59#include "proto.h"
60#include "type_parite.h"
61#include "utilitaires.h"
62
63
64
65
66 //----------------------------------------------
67 // Version Mtbl_cf
68 //----------------------------------------------
69
70namespace Lorene {
71Mtbl_cf sol_poisson_frontiere_double (const Map_af& mapping,
72 const Mtbl_cf& source, const Mtbl_cf& lim_func, const Mtbl_cf& lim_der,
73 int num_zone)
74
75{
76
77 // Verifications d'usage sur les zones
78 int nz = source.get_mg()->get_nzone() ;
79 assert (nz>1) ;
80 assert ((num_zone>0) && (num_zone<nz-1)) ;
81 assert(source.get_mg()->get_type_r(num_zone) == FIN) ;
82
83 assert (lim_func.get_mg() == source.get_mg()->get_angu()) ;
84 assert (lim_der.get_mg() == source.get_mg()->get_angu()) ;
85 assert (source.get_etat() != ETATNONDEF) ;
86 assert (lim_func.get_etat() != ETATNONDEF) ;
87 assert (lim_der.get_etat() != ETATNONDEF) ;
88
89 // Bases spectrales
90 const Base_val& base = source.base ;
91
92 // donnees sur la zone
93 int nr = source.get_mg()->get_nr(num_zone) ;
94 int nt = source.get_mg()->get_nt(num_zone) ;
95 int np = source.get_mg()->get_np(num_zone) ;;
96 int base_r ;
97 int l_quant, m_quant;
98
99 double alpha = mapping.get_alpha()[num_zone] ;
100 double beta = mapping.get_beta()[num_zone] ;
101 double echelle = beta/alpha ;
102 double facteur ;
103
104 //Rangement des valeurs intermediaires
105 Tbl *so ;
106 Tbl *sol_hom ;
107 Tbl *sol_part ;
108 Matrice *operateur ;
109 Matrice *nondege ;
110
111
112 Mtbl_cf resultat(source.get_mg(), base) ;
113 resultat.annule_hard() ;
114
115 for (int k=0 ; k<np+1 ; k++)
116 for (int j=0 ; j<nt ; j++)
117 if (nullite_plm(j, nt, k, np, base) == 1)
118 {
119 // calcul des nombres quantiques :
120 donne_lm(nz, num_zone, j, k, base, m_quant, l_quant, base_r) ;
121
122 // Construction de l'operateur
123 operateur = new Matrice(laplacien_mat
124 (nr, l_quant, echelle, 0, base_r)) ;
125
126 (*operateur) = combinaison(*operateur, l_quant, echelle, 0,
127 base_r) ;
128
129 // Operateur inversible
130 nondege = new Matrice(prepa_nondege(*operateur, l_quant,
131 echelle, 0, base_r)) ;
132
133 // Calcul DES DEUX SH
134 sol_hom = new Tbl(solh(nr, l_quant, echelle, base_r)) ;
135
136 // Calcul de la SP
137 so = new Tbl(nr) ;
138 so->set_etat_qcq() ;
139 for (int i=0 ; i<nr ; i++)
140 so->set(i) = source(num_zone, k, j, i) ;
141
142 sol_part = new Tbl (solp(*operateur, *nondege, alpha,
143 beta, *so, 0, base_r)) ;
144
145 //-------------------------------------------
146 // On est parti pour imposer la boundary
147 //-------------------------------------------
148 // Conditions de raccord type Dirichlet :
149 // Pour la sp :
150 double somme = 0 ;
151 for (int i=0 ; i<nr ; i++)
152 if (i%2 == 0)
153 somme += (*sol_part)(i) ;
154 else
155 somme -= (*sol_part)(i) ;
156
157 facteur = (lim_func(num_zone-1, k, j, 0)-somme)
158 * pow(echelle-1, l_quant+1) ;
159
160 for (int i=0 ; i<nr ; i++)
161 sol_part->set(i) +=
162 facteur*(*sol_hom)(1, i) ;
163
164 // pour l'autre solution homogene :
165 facteur = - pow(echelle-1, 2*l_quant+1) ;
166 for (int i=0 ; i<nr ; i++)
167 sol_hom->set(0, i) +=
168 facteur*(*sol_hom)(1, i) ;
169
170 // Condition de raccord de type Neumann :
171 double val_der_solp = 0 ;
172 for (int i=0 ; i<nr ; i++)
173 if (i%2 == 0)
174 val_der_solp -= i*i*(*sol_part)(i)/alpha ;
175 else
176 val_der_solp += i*i*(*sol_part)(i)/alpha ;
177
178 double val_der_solh = 0 ;
179 for (int i=0 ; i<nr ; i++)
180 if (i%2 == 0)
181 val_der_solh -= i*i*(*sol_hom)(0, i)/alpha ;
182 else
183 val_der_solh += i*i*(*sol_hom)(0, i)/alpha ;
184
185 assert (val_der_solh != 0) ;
186
187 facteur = (lim_der(num_zone-1, k, j, 0)-val_der_solp) /
188 val_der_solh ;
189
190 for (int i=0 ; i<nr ; i++)
191 sol_part->set(i) +=
192 facteur*(*sol_hom)(0, i) ;
193
194 // solp contient le bon truc (normalement ...)
195 for (int i=0 ; i<nr ; i++)
196 resultat.set(num_zone, k, j, i) = (*sol_part)(i) ;
197
198 delete operateur ;
199 delete nondege ;
200 delete so ;
201 delete sol_hom ;
202 delete sol_part ;
203 }
204 return resultat ;
205}
206}
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
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 pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene prototypes.
Definition app_hor.h:64