LORENE
ope_sec_order_mat.C
1/*
2 * Copyright (c) 2003 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 version 2
8 * as published by the Free Software Foundation.
9 *
10 * LORENE is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with LORENE; if not, write to the Free Software
17 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 *
19 */
20
21char ope_sec_order_mat_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_sec_order/ope_sec_order_mat.C,v 1.3 2014/10/13 08:53:35 j_novak Exp $" ;
22
23/*
24 * $Id: ope_sec_order_mat.C,v 1.3 2014/10/13 08:53:35 j_novak Exp $
25 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_sec_order/ope_sec_order_mat.C,v 1.3 2014/10/13 08:53:35 j_novak Exp $
26 *
27 */
28#include <cmath>
29#include <cstdlib>
30
31#include "proto.h"
32#include "ope_elementary.h"
33
34
35 //-----------------------------------
36 // Routine pour les cas non prevus --
37 //-----------------------------------
38
39namespace Lorene {
40Matrice _sec_order_mat_pas_prevu(int, double, double, double, double) {
41 cout << "Sec_order : base not implemented..." << endl ;
42 abort() ;
43 exit(-1) ;
44 Matrice res(1, 1) ;
45 return res;
46}
47
48
49
50 //-------------------------
51 //-- CAS R_CHEB -----
52 //------------------------
53
54Matrice _sec_order_mat_r_cheb (int n, double alpha,
55 double a, double b, double c) {
56
57 double* vect = new double[n] ;
58
59 Matrice dd (n,n) ;
60 dd.set_etat_qcq() ;
61 Matrice df (n,n) ;
62 df.set_etat_qcq() ;
63 Matrice ff (n,n) ;
64 ff.set_etat_qcq() ;
65
66
67 // Calcul terme en d^2...
68 for (int i=0 ; i<n ; i++) {
69 for (int j=0 ; j<n ; j++)
70 vect[j] = 0 ;
71 vect[i] = 1 ;
72
73 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
74
75 for (int j=0 ; j<n ; j++)
76 dd.set(j,i) = vect[j] ;
77 }
78
79 // Calcul terme en d...
80 for (int i=0 ; i<n ; i++) {
81 for (int j=0 ; j<n ; j++)
82 vect[j] = 0 ;
83 vect[i] = 1 ;
84
85 sxdsdx_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
86
87 for (int j=0 ; j<n ; j++)
88 df.set(j,i) = vect[j] ;
89 }
90
91 // Calcul terme Id
92 for (int i=0 ; i<n ; i++) {
93 for (int j=0 ; j<n ; j++)
94 ff.set(j,i) = 0 ;
95 ff.set(i,i) = 1 ;
96 }
97
98 delete [] vect ;
99
100 return a*dd/alpha/alpha+b*df/alpha+c*ff ;
101}
102
104 if (ope_mat != 0x0)
105 delete ope_mat ;
106
107 // Routines de derivation
110 static int nap = 0 ;
111
112 // Premier appel
113 if (nap==0) {
114 nap = 1 ;
115 for (int i=0 ; i<MAX_BASE ; i++) {
116 sec_order_mat[i] = _sec_order_mat_pas_prevu ;
117 }
118 // Les routines existantes
119 sec_order_mat[R_CHEB >> TRA_R] = _sec_order_mat_r_cheb ;
120 }
123}
124}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Matrix handling.
Definition matrice.h:152
Matrice * ope_mat
Pointer on the matrix representation of the operator.
double alpha
Parameter of the associated mapping.
int base_r
Radial basis of decomposition.
int nr
Number of radial points.
double b_param
The parameter .
virtual void do_ope_mat() const
Computes the matrix of the operator.
double a_param
The parameter .
double c_param
The parameter .
#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)
Lorene prototypes.
Definition app_hor.h:64