LORENE
multx2_1d.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 multx2_1d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/multx2_1d.C,v 1.5 2015/03/05 08:49:32 j_novak Exp $" ;
24 
25 /*
26  * $Id: multx2_1d.C,v 1.5 2015/03/05 08:49:32 j_novak Exp $
27  * $Log: multx2_1d.C,v $
28  * Revision 1.5 2015/03/05 08:49:32 j_novak
29  * Implemented operators with Legendre bases.
30  *
31  * Revision 1.4 2014/10/13 08:53:24 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.3 2014/10/06 15:16:06 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.2 2002/10/16 14:36:58 j_novak
38  * Reorganization of #include instructions of standard C++, in order to
39  * use experimental version 3 of gcc.
40  *
41  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
42  * LORENE
43  *
44  * Revision 2.0 1999/10/11 09:53:08 phil
45  * *** empty log message ***
46  *
47  *
48  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/multx2_1d.C,v 1.5 2015/03/05 08:49:32 j_novak Exp $
49  *
50  */
51 
52 // Includes
53 #include <cstdlib>
54 #include <cassert>
55 
56 #include "headcpp.h"
57 #include "type_parite.h"
58 #include "proto.h"
59 
60  //-----------------------------------
61  // Routine pour les cas non prevus --
62  //-----------------------------------
63 
64 namespace Lorene {
65 void _multx2_1d_pas_prevu(int nr, double* tb, double *res) {
66  cout << "multx2 pas prevu..." << tb << " " << res << endl ;
67  cout << "nr : " << nr << endl ;
68  abort() ;
69  exit(-1) ;
70 }
71 
72  //---------------
73  // cas R_CHEBP --
74  //---------------
75 
76 void _multx2_1d_r_chebp(int nr, double* tb, double *xo) {
77 
78  assert (nr>2) ;
79 
80  xo[0] = (2*tb[0]+tb[1])/4 ;
81  xo[1] = (2*tb[0]+2*tb[1]+tb[2])/4 ;
82 
83  for (int i=2 ; i<nr-1 ; i++)
84  xo[i] = (tb[i-1]+2*tb[i]+tb[i+1])/4 ;
85 
86  xo[nr-1] = (tb[nr-2]+2*tb[nr-1])/4 ;
87 }
88 
89 
90  //---------------
91  // cas R_CHEBI --
92  //---------------
93 
94 void _multx2_1d_r_chebi(int nr, double* tb, double *xo){
95 
96  assert(nr>2) ;
97 
98  xo[0] = (3*tb[0]+tb[1])/4 ;
99 
100  for (int i=1 ; i<nr-1 ; i++)
101  xo[i] = (tb[i-1]+2*tb[i]+tb[i+1])/4 ;
102 
103  xo[nr-1] = 0. ;//(tb[nr-2]+2*tb[nr-1])/4 ;
104 }
105 
106  //---------------
107  // cas R_CHEB --
108  //---------------
109 
110 void _multx2_1d_r_cheb(int nr, double* tb, double *xo){
111 
112  assert(nr>3) ;
113 
114  xo[0] = (2*tb[0]+tb[2])/4 ;
115  xo[1] = (3*tb[1]+tb[3])/4 ;
116  xo[2] = (2*tb[0]+2*tb[2]+tb[4])/4 ;
117 
118  for (int i=3 ; i<nr-2 ; i++)
119  xo[i] = (tb[i-2]+2*tb[i]+tb[i+2])/4 ;
120 
121  for (int i=nr-2 ; i<nr ; i++)
122  xo[i] = (tb[i-2]+2*tb[i])/4 ;
123 }
124 
125 
126  //--------------
127  // cas R_LEGP --
128  //--------------
129 
130 void _multx2_1d_r_legp(int nr, double* tb, double *xo) {
131 
132  assert (nr>2) ;
133 
134  xo[0] = tb[0]/3. + 2.*tb[1]/15. ;
135 
136  for (int i=1 ; i<nr-1 ; i++)
137  xo[i] = double(2*i*(2*i-1))*tb[i-1]/double((4*i-1)*(4*i-3))
138  + (double(4*i*i)/double((4*i-1)*(4*i+1))
139  + double((2*i+1)*(2*i+1))/double((4*i+1)*(4*i+3)))*tb[i]
140  + double((2*i+1)*(2*i+2))*tb[i+1]/double((4*i+3)*(4*i+5)) ;
141  xo[nr-1] = double((2*nr-2)*(2*nr-3))*tb[nr-2]/double((4*nr-5)*(4*nr-7))
142  + (double(4*(nr-1)*(nr-1))/double((4*nr-5)*(4*nr-3))
143  + double((2*nr-1)*(2*nr-1))/double((4*nr-3)*(4*nr-1)))*tb[nr-1] ;
144 }
145 
146 
147  //--------------
148  // cas R_LEGI --
149  //--------------
150 
151 void _multx2_1d_r_legi(int nr, double* tb, double *xo){
152 
153  assert(nr>2) ;
154 
155  xo[0] = 0.6*tb[0] + 6.*tb[1]/35. ;
156 
157  for (int i=1 ; i<nr-1 ; i++)
158  xo[i] = double(2*i*(2*i+1))*tb[i-1]/double((4*i+1)*(4*i-1))
159  + (double((2*i+1)*(2*i+1))/double((4*i+1)*(4*i+3))
160  + double((2*i+2)*(2*i+2))/double((4*i+3)*(4*i+5)) )*tb[i]
161  + double((2*i+2)*(2*i+3))*tb[i+1]/double((4*i+7)*(4*i+5)) ;
162  xo[nr-1] = 0. ;
163 }
164 
165  //-------------
166  // cas R_LEG --
167  //-------------
168 
169 void _multx2_1d_r_leg(int nr, double* tb, double *xo){
170 
171  assert(nr>3) ;
172 
173  xo[0] = tb[0]/3. + 2.*tb[2]/15. ;
174  xo[1] = 0.6*tb[1] + 6.*tb[3]/35. ;
175  for (int i=2 ; i<nr-2 ; i++)
176  xo[i] = double(i*(i-1))*tb[i-2]/double((2*i-1)*(2*i-3))
177  + (double(i*i)/double((2*i-1)*(2*i+1))
178  + double((i+1)*(i+1))/double((2*i+1)*(2*i+3)))*tb[i]
179  + double((i+1)*(i+2))*tb[i+2]/double((2*i+3)*(2*i+5)) ;
180  for (int i=nr-2 ; i<nr ; i++)
181  xo[i] = double(i*(i-1))*tb[i-2]/double((2*i-1)*(2*i-3))
182  + (double(i*i)/double((2*i-1)*(2*i+1))
183  + double((i+1)*(i+1))/double((2*i+1)*(2*i+3)))*tb[i] ;
184 }
185 
186 
187 
188 
189  //----------------------
190  // La routine a appeler
191  //----------------------
192 
193 void multx2_1d(int nr, double **tb, int base_r) // Version appliquee a this
194 {
195 
196 // Routines de derivation
197 static void (*multx2_1d[MAX_BASE])(int, double *, double*) ;
198 static int nap = 0 ;
199 
200  // Premier appel
201  if (nap==0) {
202  nap = 1 ;
203  for (int i=0 ; i<MAX_BASE ; i++) {
204  multx2_1d[i] = _multx2_1d_pas_prevu ;
205  }
206  // Les routines existantes
207  multx2_1d[R_CHEB >> TRA_R] = _multx2_1d_r_cheb ;
208  multx2_1d[R_CHEBP >> TRA_R] = _multx2_1d_r_chebp ;
209  multx2_1d[R_CHEBI >> TRA_R] = _multx2_1d_r_chebi ;
210  multx2_1d[R_LEG >> TRA_R] = _multx2_1d_r_leg ;
211  multx2_1d[R_LEGP >> TRA_R] = _multx2_1d_r_legp ;
212  multx2_1d[R_LEGI >> TRA_R] = _multx2_1d_r_legi ;
213  }
214 
215  double *result = new double[nr] ;
216  multx2_1d[base_r](nr, *tb, result) ;
217 
218  delete [] (*tb) ;
219  (*tb) = result ;
220 }
221 }
Lorene prototypes.
Definition: app_hor.h:64
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define R_LEG
base de Legendre ordinaire (fin)
Definition: type_parite.h:182
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166