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 char 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 
70 namespace Lorene {
71 Mtbl_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