LORENE
solp_helmholtz_plus.C
1 /*
2  * Copyright (c) 1999-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 solp_helmholtz_plus_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/solp_helmholtz_plus.C,v 1.5 2014/10/13 08:53:31 j_novak Exp $" ;
24 
25 /*
26  * $Id: solp_helmholtz_plus.C,v 1.5 2014/10/13 08:53:31 j_novak Exp $
27  * $Log: solp_helmholtz_plus.C,v $
28  * Revision 1.5 2014/10/13 08:53:31 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.4 2014/10/06 15:16:10 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.3 2008/02/18 13:53:45 j_novak
35  * Removal of special indentation instructions.
36  *
37  * Revision 1.2 2004/01/15 09:15:37 p_grandclement
38  * Modification and addition of the Helmholtz operators
39  *
40  * Revision 1.1 2003/12/11 14:48:49 p_grandclement
41  * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
42  *
43  *
44  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/solp_helmholtz_plus.C,v 1.5 2014/10/13 08:53:31 j_novak Exp $
45  *
46  */
47 
48 //fichiers includes
49 #include <cstdio>
50 #include <cstdlib>
51 #include <cmath>
52 
53 #include "matrice.h"
54 #include "type_parite.h"
55 #include "proto.h"
56 
57  //------------------------------------
58  // Routine pour les cas non prevus --
59  //------------------------------------
60 namespace Lorene {
61 Tbl _solp_helmholtz_plus_pas_prevu (const Matrice &, const Matrice &,
62  const Tbl &, double, double) {
63  cout << " Solution homogene pas prevue ..... : "<< endl ;
64  abort() ;
65  exit(-1) ;
66  Tbl res(1) ;
67  return res;
68 }
69 
70 
71 
72  //-------------------
73  //-- R_CHEB -----
74  //-------------------
75 Tbl _solp_helmholtz_plus_r_cheb (const Matrice &lap, const Matrice &nondege,
76  const Tbl &source, double alpha, double beta) {
77 
78  int n = lap.get_dim(0) ;
79  int dege = n-nondege.get_dim(0) ;
80  assert (dege ==2) ;
81 
82  Tbl source_aux(source*alpha*alpha) ;
83  Tbl xso(source_aux) ;
84  Tbl xxso(source_aux) ;
85  multx_1d(n, &xso.t, R_CHEB) ;
86  multx_1d(n, &xxso.t, R_CHEB) ;
87  multx_1d(n, &xxso.t, R_CHEB) ;
88  source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
89 
90  source_aux = cl_helmholtz_plus (source_aux, R_CHEB) ;
91 
92  Tbl so(n-dege) ;
93  so.set_etat_qcq() ;
94  for (int i=0 ; i<n-dege ; i++)
95  so.set(i) = source_aux(i) ;
96 
97  Tbl auxi(nondege.inverse(so)) ;
98 
99  Tbl res(n) ;
100  res.set_etat_qcq() ;
101  for (int i=dege ; i<n ; i++)
102  res.set(i) = auxi(i-dege) ;
103 
104  for (int i=0 ; i<dege ; i++)
105  res.set(i) = 0 ;
106  return res ;
107 }
108 
109  //-------------------
110  //-- R_CHEBP ----
111  //-------------------
112 Tbl _solp_helmholtz_plus_r_chebp (const Matrice &lap, const Matrice &nondege,
113  const Tbl &source, double alpha, double) {
114 
115  int n = lap.get_dim(0) ;
116  int dege = n-nondege.get_dim(0) ;
117  assert (dege ==1) ;
118 
119  Tbl source_aux(source*alpha*alpha) ;
120  source_aux = cl_helmholtz_plus (source_aux, R_CHEBP) ;
121 
122  Tbl so(n-dege) ;
123  so.set_etat_qcq() ;
124  for (int i=0 ; i<n-dege ; i++)
125  so.set(i) = source_aux(i) ;
126 
127  Tbl auxi(nondege.inverse(so)) ;
128 
129  Tbl res(n) ;
130  res.set_etat_qcq() ;
131  for (int i=dege ; i<n ; i++)
132  res.set(i) = auxi(i-dege) ;
133 
134  for (int i=0 ; i<dege ; i++)
135  res.set(i) = 0 ;
136  return res ;
137 }
138  //-------------------
139  //-- Fonction ----
140  //-------------------
141 
142 
143 Tbl solp_helmholtz_plus (const Matrice &lap, const Matrice &nondege,
144  const Tbl &source, double alpha, double beta,
145  int base_r) {
146 
147  // Routines de derivation
148  static Tbl (*solp_helmholtz_plus[MAX_BASE]) (const Matrice&, const Matrice&,
149  const Tbl&, double, double) ;
150  static int nap = 0 ;
151 
152  // Premier appel
153  if (nap==0) {
154  nap = 1 ;
155  for (int i=0 ; i<MAX_BASE ; i++) {
156  solp_helmholtz_plus[i] = _solp_helmholtz_plus_pas_prevu ;
157  }
158  // Les routines existantes
159  solp_helmholtz_plus[R_CHEB >> TRA_R] = _solp_helmholtz_plus_r_cheb ;
160  solp_helmholtz_plus[R_CHEBP >> TRA_R] = _solp_helmholtz_plus_r_chebp ;
161  }
162 
163  Tbl res(solp_helmholtz_plus[base_r] (lap, nondege, source, alpha, beta)) ;
164  return res ;
165 }
166 }
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:403
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Lorene prototypes.
Definition: app_hor.h:64