LORENE
comb_lin_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
23char comb_lin_helmholtz_plusC[] = "$Header $" ;
24
25/*
26 * $Id: comb_lin_helmholtz_plus.C,v 1.5 2014/10/13 08:53:28 j_novak Exp $
27 * $Log: comb_lin_helmholtz_plus.C,v $
28 * Revision 1.5 2014/10/13 08:53:28 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.4 2014/10/06 15:16:08 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.3 2008/02/18 13:53:43 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/comb_lin_helmholtz_plus.C,v 1.5 2014/10/13 08:53:28 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// Version Matrice --> Matrice
59namespace Lorene {
60Matrice _cl_helmholtz_plus_pas_prevu (const Matrice& so) {
61 cout << "CL Helmholtz plus not implemented" << endl ;
62 abort() ;
63 exit(-1) ;
64 return so;
65}
66
67
68 //-------------------
69 //-- R_CHEBP -----
70 //-------------------
71
72
73Matrice _cl_helmholtz_plus_r_chebp (const Matrice &source) {
74
75 int n = source.get_dim(0) ;
76 assert (n == source.get_dim(1)) ;
77
78 Matrice barre(source) ;
79
80 int dirac = 1 ;
81 for (int i=0 ; i<n-2 ; i++) {
82 for (int j=0 ; j<n ; j++)
83 barre.set(i, j) = (1+dirac)*source(i, j)-source(i+2, j) ;
84 if (i==0) dirac = 0 ;
85 }
86
87 Matrice tilde(barre) ;
88 for (int i=0 ; i<n-4 ; i++)
89 for (int j=0 ; j<n ; j++)
90 tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
91
92 Matrice res(tilde) ;
93 for (int i=0 ; i<n-4 ; i++)
94 for (int j=0 ; j<n ; j++)
95 res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
96
97 return res ;
98}
99
100
101 //-------------------
102 //-- R_CHEB ------
103 //-------------------
104
105Matrice _cl_helmholtz_plus_r_cheb (const Matrice& source) {
106
107
108 int n = source.get_dim(0) ;
109 assert (n==source.get_dim(1)) ;
110
111 Matrice barre(source) ;
112 int dirac = 1 ;
113 for (int i=0 ; i<n-2 ; i++) {
114 for (int j=0 ; j<n ; j++)
115 barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j))
116 /(i+1) ;
117 if (i==0) dirac = 0 ;
118 }
119
120 Matrice res(barre) ;
121 for (int i=0 ; i<n-4 ; i++)
122 for (int j=0 ; j<n ; j++)
123 res.set(i, j) = barre(i, j)-barre(i+2, j) ;
124
125 return res ;
126}
127
128
129 //-------------------------
130 //- La routine a appeler ---
131 //---------------------------
132
133Matrice cl_helmholtz_plus (const Matrice &source, int base_r) {
134
135 // Routines de derivation
136 static Matrice (*cl_helmholtz_plus[MAX_BASE]) (const Matrice &) ;
137 static int nap = 0 ;
138
139 // Premier appel
140 if (nap==0) {
141 nap = 1 ;
142 for (int i=0 ; i<MAX_BASE ; i++) {
143 cl_helmholtz_plus[i] = _cl_helmholtz_plus_pas_prevu ;
144 }
145 // Les routines existantes
146 cl_helmholtz_plus[R_CHEB >> TRA_R] = _cl_helmholtz_plus_r_cheb ;
147 cl_helmholtz_plus[R_CHEBP >> TRA_R] = _cl_helmholtz_plus_r_chebp ;
148 }
149
150 Matrice res(cl_helmholtz_plus[base_r](source)) ;
151 return res ;
152}
153
154
155//************************ TBL Versions *************************************
156
157
158
159
160Tbl _cl_helmholtz_plus_pas_prevu (const Tbl &so) {
161
162 cout << "Linear combination for Helmholtz plus not implemented..." << endl ;
163 abort() ;
164 exit(-1) ;
165 return so;
166}
167
168
169 //-------------------
170 //-- R_CHEBP -------
171 //--------------------
172Tbl _cl_helmholtz_plus_r_chebp (const Tbl& source) {
173
174 int n = source.get_dim(0) ;
175
176 Tbl barre(source) ;
177 int dirac = 1 ;
178 for (int i=0 ; i<n-2 ; i++) {
179 barre.set(i) = (1+dirac)*source(i)-source(i+2) ;
180 if (i==0) dirac = 0 ;
181 }
182
183 Tbl tilde(barre) ;
184 for (int i=0 ; i<n-4 ; i++)
185 tilde.set(i) = barre(i)-barre(i+2) ;
186
187 Tbl res(tilde) ;
188 for (int i=0 ; i<n-4 ; i++)
189 res.set(i) = tilde(i)-tilde(i+1) ;
190
191 return res ;
192}
193
194 //-------------------
195 //-- R_CHEB -------
196 //--------------------
197Tbl _cl_helmholtz_plus_r_cheb (const Tbl& source) {
198
199 int n = source.get_dim(0) ;
200
201 Tbl barre(source) ;
202 int dirac = 1 ;
203 for (int i=0 ; i<n-2 ; i++) {
204 barre.set(i) = ((1+dirac)*source(i)-source(i+2))
205 /(i+1) ;
206 if (i==0) dirac = 0 ;
207 }
208
209 Tbl res(barre) ;
210 for (int i=0 ; i<n-4 ; i++)
211 res.set(i) = barre(i)-barre(i+2) ;
212
213 return res ;
214}
215 //----------------------------
216 //- Routine a appeler ---
217 //------------------------------
218
219Tbl cl_helmholtz_plus (const Tbl &source, int base_r) {
220
221 // Routines de derivation
222 static Tbl (*cl_helmholtz_plus[MAX_BASE])(const Tbl &) ;
223 static int nap = 0 ;
224
225 // Premier appel
226 if (nap==0) {
227 nap = 1 ;
228 for (int i=0 ; i<MAX_BASE ; i++) {
229 cl_helmholtz_plus[i] = _cl_helmholtz_plus_pas_prevu ;
230 }
231 // Les routines existantes
232 cl_helmholtz_plus[R_CHEB >> TRA_R] = _cl_helmholtz_plus_r_cheb ;
233 cl_helmholtz_plus[R_CHEBP >> TRA_R] = _cl_helmholtz_plus_r_chebp ;
234 }
235
236 Tbl res(cl_helmholtz_plus[base_r](source)) ;
237 return res ;
238}
239}
int get_dim(int i) const
Returns the dimension of the matrix.
Definition matrice.C:260
#define MAX_BASE
Nombre max. de bases differentes.
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64