LORENE
raccord_c1.C
1/*
2 * Copyright (c) 2000-2001 Eric Gourgoulhon
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 raccord_c1_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/raccord_c1.C,v 1.4 2014/10/13 08:53:32 j_novak Exp $" ;
24
25/*
26 * $Id: raccord_c1.C,v 1.4 2014/10/13 08:53:32 j_novak Exp $
27 * $Log: raccord_c1.C,v $
28 * Revision 1.4 2014/10/13 08:53:32 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.3 2010/01/26 16:47:40 e_gourgoulhon
32 * Added the Scalar version.
33 *
34 * Revision 1.2 2003/12/19 16:21:49 j_novak
35 * Shadow hunt
36 *
37 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
38 * LORENE
39 *
40 * Revision 2.0 2000/03/22 10:08:55 eric
41 * *** empty log message ***
42 *
43 *
44 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/raccord_c1.C,v 1.4 2014/10/13 08:53:32 j_novak Exp $
45 *
46 */
47
48// Headers Lorene
49#include "cmp.h"
50#include "scalar.h"
51
52namespace Lorene {
53Cmp raccord_c1(const Cmp& uu, int l1) {
54
55 const Map_radial* mpi = dynamic_cast<const Map_radial*>( uu.get_mp() ) ;
56
57 if (mpi == 0x0) {
58 cout <<
59 "raccord_c1 : The mapping does not belong to the class Map_radial !"
60 << endl ;
61 abort() ;
62 }
63
64 const Map_radial& mp = *mpi ;
65
66 const Mg3d& mg = *(mp.get_mg()) ;
67
68#ifndef NDEBUG
69 int nz = mg.get_nzone() ;
70#endif
71 assert(l1 > 0) ;
72 assert(l1 < nz-1) ;
73
74 int l0 = l1 - 1 ; // index of left domain
75 int l2 = l1 + 1 ; // index of right domain
76
77
78 Mtbl dxdr = mp.dxdr ;
79 Mtbl r2 = mp.r * mp.r ;
80
81 Cmp resu = uu ;
82
83 Valeur& va = resu.va ;
84
85 va.coef_i() ;
86 va.set_etat_c_qcq() ;
87
88 va.c->set_etat_qcq() ;
89 va.c->t[l1]->set_etat_qcq() ;
90
91 int np = mg.get_np(l1) ;
92 int nt = mg.get_nt(l1) ;
93 assert( mg.get_np(l0) == np ) ;
94 assert( mg.get_nt(l0) == nt ) ;
95 assert( mg.get_np(l2) == np ) ;
96 assert( mg.get_nt(l2) == nt ) ;
97
98 int nr0 = mg.get_nr(l0) ;
99 int nr1 = mg.get_nr(l1) ;
100
101 for (int k=0; k<np; k++) {
102 for (int j=0; j<nt; j++) {
103 double lam0 = uu(l0, k, j, nr0-1) ;
104 double lam1 = uu.dsdr()(l0, k, j, nr0-1) / dxdr(l1, k, j, 0) ;
105 double mu0 = uu(l2, k, j, 0) ;
106 double mu1 = uu.dsdr()(l2, k, j, 0) / dxdr(l1, k, j, nr1-1) ;
107
108 if ( mg.get_type_r(l2) == UNSURR ) {
109 mu1 /= r2(l2, k, j, 0) ;
110 }
111
112 double s0 = 0.25 * (mu0 + lam0) ;
113 double s1 = 0.25 * (mu1 + lam1) ;
114 double d0 = 0.25 * (mu0 - lam0) ;
115 double d1 = 0.25 * (mu1 - lam1) ;
116
117 double a0 = 2. * s0 - d1 ;
118 double a1 = 3. * d0 - s1 ;
119 double a2 = d1 ;
120 double a3 = s1 - d0 ;
121
122 for (int i=0; i<nr1; i++) {
123 double x = mg.get_grille3d(l1)->x[i] ;
124 va.set(l1, k, j, i) = a0 + x * ( a1 + x * ( a2 + x * a3 ) ) ;
125 }
126
127 }
128 }
129
130 return resu ;
131
132}
133
134/*
135 * Scalar version
136 */
137Scalar raccord_c1(const Scalar& uu, int l1) {
138
139 Cmp cuu(uu) ;
140 return Scalar( raccord_c1(cuu, l1) ) ;
141
142}
143}
Lorene prototypes.
Definition app_hor.h:64