LORENE
ope_sec_order_r2_solp.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_solp.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_solp.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  // Cl version Tbl -> Tbl --
37  //------------------------------------
38 namespace Lorene {
39 Tbl _cl_sec_order_r2_pas_prevu (const Tbl &so) {
40 
41  cout << "Linear combination for Sec_order_r2 not implemented..." << endl ;
42  abort() ;
43  exit(-1) ;
44  return so;
45 }
46 
47  //-------------------
48  //-- R_CHEB -------
49  //--------------------
50 Tbl _cl_sec_order_r2_r_cheb (const Tbl& source) {
51 
52  int n = source.get_dim(0) ;
53 
54  Tbl barre(source) ;
55  int dirac = 1 ;
56  for (int i=0 ; i<n-2 ; i++) {
57  barre.set(i) = ((1+dirac)*source(i)-source(i+2))
58  /(i+1) ;
59  if (i==0) dirac = 0 ;
60  }
61 
62  Tbl res(barre) ;
63  for (int i=0 ; i<n-4 ; i++)
64  res.set(i) = barre(i)-barre(i+2) ;
65 
66  return res ;
67 }
68 
69 
70  //----------------------------
71  //- Routine a appeler ---
72  //------------------------------
73 
74 Tbl cl_sec_order_r2 (const Tbl &source, int base_r) {
75 
76  // Routines de derivation
77  static Tbl (*cl_sec_order_r2[MAX_BASE])(const Tbl &) ;
78  static int nap = 0 ;
79 
80  // Premier appel
81  if (nap==0) {
82  nap = 1 ;
83  for (int i=0 ; i<MAX_BASE ; i++) {
84  cl_sec_order_r2[i] = _cl_sec_order_r2_pas_prevu ;
85  }
86  // Les routines existantes
87  cl_sec_order_r2[R_CHEB >> TRA_R] = _cl_sec_order_r2_r_cheb ;
88  }
89 
90  Tbl res(cl_sec_order_r2[base_r](source)) ;
91  return res ;
92 }
93 
94 
95  //*******************************
96  // CALCUL SP proprement parler
97  //*******************************
98 
99  //------------------------------------
100  // Routine pour les cas non prevus --
101  //------------------------------------
102 Tbl _solp_sec_order_r2_pas_prevu (const Matrice &, const Matrice &,
103  const Tbl &) {
104  cout << " Solution particuliere pas prevue in sec_order_r2..... : "<< endl ;
105  abort() ;
106  exit(-1) ;
107  Tbl res(1) ;
108  return res;
109 }
110 
111  //-------------------
112  //-- R_CHEB -----
113  //-------------------
114 
115 Tbl _solp_sec_order_r2_r_cheb (const Matrice &lap, const Matrice &nondege,
116  const Tbl &source) {
117 
118  int n = lap.get_dim(0) ;
119  int dege = n-nondege.get_dim(0) ;
120  assert (dege ==2) ;
121 
122  Tbl source_aux (cl_sec_order_r2 (source, R_CHEB)) ;
123 
124  Tbl so(n-dege) ;
125  so.set_etat_qcq() ;
126  for (int i=0 ; i<n-dege ; i++)
127  so.set(i) = source_aux(i) ;
128 
129  Tbl auxi(nondege.inverse(so)) ;
130 
131  Tbl res(n) ;
132  res.set_etat_qcq() ;
133  for (int i=dege ; i<n ; i++)
134  res.set(i) = auxi(i-dege) ;
135 
136  for (int i=0 ; i<dege ; i++)
137  res.set(i) = 0 ;
138  return res ;
139 }
140 
141 
142 
143 Tbl Ope_sec_order_r2::get_solp (const Tbl& so) const {
144 
145  if (non_dege == 0x0)
146  do_non_dege() ;
147 
148  // Routines de derivation
149  static Tbl (*solp_sec_order_r2[MAX_BASE]) (const Matrice&, const Matrice&,
150  const Tbl&) ;
151  static int nap = 0 ;
152 
153  // Premier appel
154  if (nap==0) {
155  nap = 1 ;
156  for (int i=0 ; i<MAX_BASE ; i++) {
157  solp_sec_order_r2[i] = _solp_sec_order_r2_pas_prevu ;
158  }
159  // Les routines existantes
160  solp_sec_order_r2[R_CHEB >> TRA_R] = _solp_sec_order_r2_r_cheb ;
161  }
162 
163  Tbl res(solp_sec_order_r2[base_r] (*ope_mat, *non_dege, so)) ;
164 
165  Tbl valeurs (val_solp (res, alpha, base_r)) ;
166  sp_plus = valeurs(0) ;
167  sp_minus = valeurs(1) ;
168  dsp_plus = valeurs(2) ;
169  dsp_minus = valeurs(3) ;
170 
171  return res ;
172 }
173 }
double alpha
Parameter of the associated mapping.
Lorene prototypes.
Definition: app_hor.h:67
Matrice * ope_mat
Pointer on the matrix representation of the operator.
double dsp_minus
Value of the derivative of the particular solution at the inner boundary.
double sp_minus
Value of the particular solution at the inner boundary.
int base_r
Radial basis of decomposition.
double sp_plus
Value of the particular solution at the outer boundary.
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
virtual Tbl get_solp(const Tbl &so) const
Computes the particular solution, given the source so .
double dsp_plus
Value of the derivative of the particular solution at the outer boundary.
Matrix handling.
Definition: matrice.h:152
virtual void do_non_dege() const
Computes the non-degenerated matrix of the operator.
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:423
int get_dim(int i) const
Returns the dimension of the matrix.
Definition: matrice.C:263
Basic array class.
Definition: tbl.h:164
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
Matrice * non_dege
Pointer on the non-degenerated matrix of the operator.
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166