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 
23 char 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 
52 namespace Lorene {
53 Cmp 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  */
137 Scalar 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