LORENE
ope_sec_order_r2_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 
21 
22 
23 /*
24  * $Id: ope_sec_order_r2_mat.C,v 1.4 2016/12/05 16:18:13 j_novak Exp $
25  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_sec_order_r2/ope_sec_order_r2_mat.C,v 1.4 2016/12/05 16:18:13 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 
39 namespace Lorene {
40 Matrice _sec_order_r2_mat_pas_prevu(int, double, double, double,
41  double, double) {
42  cout << "Sec_order_r2 : base not implemented..." << endl ;
43  abort() ;
44  exit(-1) ;
45  Matrice res(1, 1) ;
46  return res;
47 }
48 
49 
50 
51  //-------------------------
52  //-- CAS R_CHEB -----
53  //------------------------
54 
55 Matrice _sec_order_r2_mat_r_cheb (int n, double alpha, double beta,
56  double a, double b, double c) {
57 
58  double echelle = beta / alpha ;
59 
60  double* vect = new double[n] ;
61  double* auxi = new double[n] ;
62  double* res = new double[n] ;
63 
64  Matrice dd (n,n) ;
65  dd.set_etat_qcq() ;
66  Matrice df (n,n) ;
67  df.set_etat_qcq() ;
68  Matrice ff (n,n) ;
69  ff.set_etat_qcq() ;
70 
71 
72  // Calcul terme en R2 d^2...
73  for (int i=0 ; i<n ; i++) {
74  for (int j=0 ; j<n ; j++)
75  vect[j] = 0 ;
76  vect[i] = 1 ;
77 
78  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
79 
80  for (int j=0 ; j<n ; j++)
81  auxi[j] = vect[j] ;
82  multx_1d (n, &auxi, R_CHEB) ;
83  multx_1d (n, &auxi, R_CHEB) ;
84  for (int j=0 ; j<n ; j++)
85  res[j] = auxi[j] ;
86 
87  for (int j=0 ; j<n ; j++)
88  auxi[j] = vect[j] ;
89  multx_1d (n, &auxi, R_CHEB) ;
90  for (int j=0 ; j<n ; j++)
91  res[j] += auxi[j]*2*echelle ;
92 
93  for (int j=0 ; j<n ; j++)
94  res[j] += echelle*echelle*vect[j] ;
95 
96  for (int j=0 ; j<n ; j++)
97  dd.set(j,i) = res[j] ;
98  }
99 
100  // Calcul terme en R d...
101  for (int i=0 ; i<n ; i++) {
102  for (int j=0 ; j<n ; j++)
103  vect[j] = 0 ;
104  vect[i] = 1 ;
105 
106  sxdsdx_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
107 
108  for (int j=0 ; j<n ; j++)
109  auxi[j] = vect[j] ;
110  multx_1d (n, &auxi, R_CHEB) ;
111  for (int j=0 ; j<n ; j++)
112  res[j] = auxi[j] ;
113 
114  for (int j=0 ; j<n ; j++)
115  res[j] += echelle*vect[j] ;
116 
117  for (int j=0 ; j<n ; j++)
118  df.set(j,i) = res[j] ;
119  }
120 
121  // Calcul terme en R d...
122  for (int i=0 ; i<n ; i++) {
123  for (int j=0 ; j<n ; j++)
124  ff.set(j,i) = 0 ;
125  ff.set(i,i) = 1 ;
126  }
127 
128  delete [] vect ;
129  delete [] auxi ;
130  delete [] res ;
131 
132  return a*dd+b*df+c*ff ;
133 }
134 
136  if (ope_mat != 0x0)
137  delete ope_mat ;
138 
139  // Routines de derivation
140  static Matrice (*sec_order_r2_mat[MAX_BASE])(int, double, double, double,
141  double, double);
142  static int nap = 0 ;
143 
144  // Premier appel
145  if (nap==0) {
146  nap = 1 ;
147  for (int i=0 ; i<MAX_BASE ; i++) {
148  sec_order_r2_mat[i] = _sec_order_r2_mat_pas_prevu ;
149  }
150  // Les routines existantes
151  sec_order_r2_mat[R_CHEB >> TRA_R] = _sec_order_r2_mat_r_cheb ;
152  }
153  ope_mat = new Matrice(sec_order_r2_mat[base_r](nr, alpha, beta, a_param,
154  b_param, c_param)) ;
155 }
156 }
double alpha
Parameter of the associated mapping.
double beta
Parameter of the associated mapping.
Lorene prototypes.
Definition: app_hor.h:67
Matrice * ope_mat
Pointer on the matrix representation of the operator.
double a_param
The parameter a .
virtual void do_ope_mat() const
Computes the matrix of the operator.
int base_r
Radial basis of decomposition.
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
Matrix handling.
Definition: matrice.h:152
double b_param
The parameter b .
int nr
Number of radial points.
double c_param
The parameter c .
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166