LORENE
cmp_raccord_zec.C
1/*
2 * Copyright (c) 2000-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 cmp_raccord_zec_C[] = "$Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_raccord_zec.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $" ;
24
25/*
26 * $Id: cmp_raccord_zec.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $
27 * $Log: cmp_raccord_zec.C,v $
28 * Revision 1.4 2014/10/13 08:52:48 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.3 2014/10/06 15:13:04 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.2 2003/10/03 15:58:45 j_novak
35 * Cleaning of some headers
36 *
37 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
38 * LORENE
39 *
40 * Revision 2.7 2001/03/30 13:38:32 phil
41 * *** empty log message ***
42 *
43 * Revision 2.6 2001/03/22 10:25:01 phil
44 * changement complet : cas plus general
45 *
46 * Revision 2.5 2001/02/08 14:21:32 phil
47 * correction de raccord_zec.C (on prend en compte le dernier coef ...)
48 *
49 * Revision 2.4 2001/01/02 11:25:37 phil
50 * *** empty log message ***
51 *
52 * Revision 2.3 2000/12/13 14:59:18 phil
53 * *** empty log message ***
54 *
55 * Revision 2.2 2000/12/13 14:49:54 phil
56 * changement nom variable appel
57 * /
58 *
59 * Revision 2.1 2000/12/13 14:12:29 phil
60 * correction bugs
61 *
62 * Revision 2.0 2000/12/13 14:09:31 phil
63 * *** empty log message ***
64 *
65 *
66 * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_raccord_zec.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $
67 *
68 */
69
70//standard
71#include <cstdlib>
72#include <cmath>
73
74// LORENE
75#include "matrice.h"
76#include "cmp.h"
77#include "proto.h"
78
79// Fait le raccord C1 dans la zec ...
80namespace Lorene {
81// Suppose (pour le moment, le meme nbre de points sur les angles ...)
82// et que la zone precedente est une coquille
83
84void Cmp::raccord_c1_zec (int puis, int nbre, int lmax) {
85
86 assert (nbre>0) ;
87 assert (etat != ETATNONDEF) ;
88 if (etat == ETATZERO)
89 return ;
90
91 // Le mapping doit etre affine :
92 const Map_af* map = dynamic_cast<const Map_af*>(mp) ;
93 if (map == 0x0) {
94 cout << "Le mapping doit etre affine" << endl ;
95 abort() ;
96 }
97
98 int nz = map->get_mg()->get_nzone() ;
99 int nr = map->get_mg()->get_nr (nz-1) ;
100 int nt = map->get_mg()->get_nt (nz-1) ;
101 int np = map->get_mg()->get_np (nz-1) ;
102
103 double alpha = map->get_alpha()[nz-1] ;
104 double r_cont = -1./2./alpha ; //Rayon de debut de la zec.
105
106 // On calcul les coefficients des puissances de 1./r
107 Tbl coef (nbre+2*lmax, nr) ;
108 coef.set_etat_qcq() ;
109
110 int* deg = new int[3] ;
111 deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
112 double* auxi = new double[nr] ;
113
114 for (int conte=0 ; conte<nbre+2*lmax ; conte++) {
115 for (int i=0 ; i<nr ; i++)
116 auxi[i] = pow(-1-cos(M_PI*i/(nr-1)), (conte+puis)) ;
117
118 cfrcheb(deg, deg, auxi, deg, auxi) ;
119 for (int i=0 ; i<nr ; i++)
120 coef.set(conte, i) = auxi[i]*pow (alpha, conte+puis) ;
121 }
122
123 delete[] deg ;
124 // Maintenant on va calculer les valeurs de la ieme derivee :
125 Tbl valeurs (nbre, nt, np+1) ;
126 valeurs.set_etat_qcq() ;
127
128 Cmp courant (*this) ;
129 double* res_val = new double[1] ;
130
131 for (int conte=0 ; conte<nbre ; conte++) {
132
133 courant.va.coef() ;
134 courant.va.ylm() ;
135 courant.va.c_cf->t[nz-1]->annule_hard() ;
136
137 int base_r = courant.va.base.get_base_r(nz-2) ;
138 for (int k=0 ; k<np+1 ; k++)
139 for (int j=0 ; j<nt ; j++)
140 if (nullite_plm(j, nt, k, np, courant.va.base) == 1) {
141
142 for (int i=0 ; i<nr ; i++)
143 auxi[i] = (*courant.va.c_cf)(nz-2, k, j, i) ;
144
145 switch (base_r) {
146 case R_CHEB :
147 som_r_cheb (auxi, nr, 1, 1, 1, res_val) ;
148 break ;
149 default :
150 cout << "Cas non prevu dans raccord_zec" << endl ;
151 abort() ;
152 break ;
153 }
154 valeurs.set(conte, k, j) = res_val[0] ;
155 }
156 Cmp copie (courant) ;
157 copie.dec2_dzpuis() ;
158 courant = copie.dsdr() ;
159 }
160
161 delete [] auxi ;
162 delete [] res_val ;
163
164 // On boucle sur les harmoniques : construction de la matrice
165 // et du second membre
166 va.coef() ;
167 va.ylm() ;
168 va.c_cf->t[nz-1]->annule_hard() ;
170
171 const Base_val& base = va.base ;
172 int base_r, l_quant, m_quant ;
173 for (int k=0 ; k<np+1 ; k++)
174 for (int j=0 ; j<nt ; j++)
175 if (nullite_plm(j, nt, k, np, va.base) == 1) {
176
177 donne_lm (nz, nz-1, j, k, base, m_quant, l_quant, base_r) ;
178
179 if (l_quant<=lmax) {
180
181 Matrice systeme (nbre, nbre) ;
182 systeme.set_etat_qcq() ;
183
184 for (int col=0 ; col<nbre ; col++)
185 for (int lig=0 ; lig<nbre ; lig++) {
186
187 int facteur = (lig%2==0) ? 1 : -1 ;
188 for (int conte=0 ; conte<lig ; conte++)
189 facteur *= puis+col+conte+2*l_quant ;
190 systeme.set(lig, col) = facteur/pow(r_cont, puis+col+lig+2*l_quant) ;
191 }
192
193 systeme.set_band(nbre, nbre) ;
194 systeme.set_lu() ;
195
196 Tbl sec_membre (nbre) ;
197 sec_membre.set_etat_qcq() ;
198 for (int conte=0 ; conte<nbre ; conte++)
199 sec_membre.set(conte) = valeurs(conte, k, j) ;
200
201 Tbl inv (systeme.inverse(sec_membre)) ;
202
203 for (int conte=0 ; conte<nbre ; conte++)
204 for (int i=0 ; i<nr ; i++)
205 va.c_cf->set(nz-1, k, j, i)+=
206 inv(conte)*coef(conte+2*l_quant, i) ;
207 }
208 else for (int i=0 ; i<nr ; i++)
209 va.c_cf->set(nz-1, k, j, i)
210 = 0 ;
211 }
212
213 va.ylm_i() ;
214 set_dzpuis (0) ;
215}
216}
Bases of the spectral expansions.
Definition base_val.h:322
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
Definition base_val.h:400
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
const Map * mp
Reference mapping.
Definition cmp.h:451
void raccord_c1_zec(int puis, int nbre, int lmax)
Performs the matching of the external domain with respect to the last shell using function like wit...
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition cmp.h:454
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:654
void dec2_dzpuis()
Decreases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
const Cmp & dsdr() const
Returns of *this .
Definition cmp_deriv.C:84
Affine radial mapping.
Definition map.h:2027
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Matrix handling.
Definition matrice.h:152
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition matrice.C:175
double & set(int j, int i)
Read/write of a particuliar element.
Definition matrice.h:277
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition matrice.C:424
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition matrice.C:364
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition matrice.C:392
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition mtbl_cf.h:205
Basic array class.
Definition tbl.h:161
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:372
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition valeur.C:712
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
void ylm_i()
Inverse of ylm()
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
#define R_CHEB
base de Chebychev ordinaire (fin)
Lorene prototypes.
Definition app_hor.h:64