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