LORENE
regularise_shift.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 
23 char regularise_shift_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/regularise_shift.C,v 1.7 2014/10/13 08:52:41 j_novak Exp $" ;
24 
25 /*
26  * $Id: regularise_shift.C,v 1.7 2014/10/13 08:52:41 j_novak Exp $
27  * $Log: regularise_shift.C,v $
28  * Revision 1.7 2014/10/13 08:52:41 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.6 2014/10/06 15:12:59 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.5 2006/04/27 09:12:32 p_grandclement
35  * First try at irrotational black holes
36  *
37  * Revision 1.4 2005/08/29 15:10:14 p_grandclement
38  * Addition of things needed :
39  * 1) For BBH with different masses
40  * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
41  * WORKING YET !!!
42  *
43  * Revision 1.3 2003/10/03 15:58:44 j_novak
44  * Cleaning of some headers
45  *
46  * Revision 1.2 2003/02/13 16:40:25 p_grandclement
47  * Addition of various things for the Bin_ns_bh project, non of them being
48  * completely tested
49  *
50  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
51  * LORENE
52  *
53  * Revision 2.6 2001/05/07 09:11:48 phil
54  * *** empty log message ***
55  *
56  * Revision 2.5 2001/01/29 14:30:48 phil
57  * ajout type rotation
58  *
59  * Revision 2.4 2000/11/02 10:18:05 phil
60  * modification du degre de l;a fonction\
61  * de correction. On passe de 2 a 3
62  *
63  * Revision 2.3 2000/10/31 10:04:28 phil
64  * correction importation
65  *
66  * Revision 2.2 2000/10/26 12:34:17 phil
67  * *** empty log message ***
68  *
69  * Revision 2.1 2000/10/26 12:31:55 phil
70  * correction orientation pour calcul du shift total
71  *
72  * Revision 2.0 2000/10/19 10:08:03 phil
73  * *** empty log message ***
74  *
75  *
76  * $Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/regularise_shift.C,v 1.7 2014/10/13 08:52:41 j_novak Exp $
77  *
78  */
79 
80 
81 //Standard
82 #include <cstdlib>
83 #include <cmath>
84 
85 //Lorene
86 #include "nbr_spx.h"
87 #include "tenseur.h"
88 
89 namespace Lorene {
90 double regle (Tenseur& shift_auto, const Tenseur& shift_comp, double omega, double omega_local) {
91 
92  Tenseur shift_old (shift_auto) ;
93 
94  double orientation = shift_auto.get_mp()->get_rot_phi() ;
95  assert ((orientation==0) || (orientation == M_PI)) ;
96  double orientation_autre = shift_comp.get_mp()->get_rot_phi() ;
97  assert ((orientation_autre==0) || (orientation_autre == M_PI)) ;
98 
99  int alignes = (orientation == orientation_autre) ? 1 : -1 ;
100 
101  // Cas triades identiques
102  if (*shift_comp.get_triad() == *shift_auto.get_triad())
103  alignes = 1 ;
104 
105  int nz = shift_auto.get_mp()->get_mg()->get_nzone() ;
106  int np = shift_auto.get_mp()->get_mg()->get_np(1) ;
107  int nt = shift_auto.get_mp()->get_mg()->get_nt(1) ;
108  int nr = shift_auto.get_mp()->get_mg()->get_nr(1) ;
109 
110  // On minimise la valeur de la derivee de B sur R :
111  Tenseur shift_tot (*shift_auto.get_mp(), 1, CON, *shift_auto.get_triad()) ;
112  shift_tot.set_etat_qcq() ;
113  shift_tot.set(0).import_asymy (alignes*shift_comp(0)) ;
114  shift_tot.set(1).import_symy (alignes*shift_comp(1)) ;
115  shift_tot.set(2).import_asymy (shift_comp(2)) ;
116 
117  shift_tot = shift_tot + shift_auto ;
118 
119  double indic = (orientation == 0) ? 1 : -1 ;
120 
121  Mtbl xa_mtbl (shift_tot.get_mp()->get_mg()) ;
122  xa_mtbl = shift_tot.get_mp()->xa ;
123  Mtbl ya_mtbl (shift_tot.get_mp()->get_mg()) ;
124  ya_mtbl = shift_tot.get_mp()->ya ;
125 
126  Tenseur tbi (shift_tot) ;
127  if (omega != 0) {
128  for (int i=0 ; i<3 ; i++) {
129  tbi.set(i).va.coef_i() ;
130  tbi.set(i).va.set_etat_c_qcq() ;
131  }
132 
133  tbi.set(0) = *shift_tot(0).va.c - indic *omega * ya_mtbl(0,0,0,0) - indic*omega_local* shift_tot.get_mp()->y ;
134  tbi.set(1) = *shift_tot(1).va.c + indic *omega * xa_mtbl(0,0,0,0) + indic*omega_local* shift_tot.get_mp()->x ;
135  tbi.set_std_base() ;
136  tbi.set(0).annule(nz-1) ;
137  tbi.set(1).annule(nz-1) ;
138  }
139 
140  Tenseur derive_r (*shift_auto.get_mp(), 1, CON, *shift_auto.get_triad()) ;
141  derive_r.set_etat_qcq() ;
142  for (int i=0 ; i<3 ; i++) {
143  if (tbi(i).get_etat() != ETATZERO)
144  derive_r.set(i) = tbi(i).dsdr() ;
145  else
146  derive_r.set(i).set_etat_zero() ;
147  }
148 
149  // On enleve un fonction pour rendre Kij regulier !
150 
151  Valeur val_hor (shift_auto.get_mp()->get_mg()) ;
152  Valeur fonction_radiale (shift_auto.get_mp()->get_mg()) ;
153  Cmp enleve (shift_auto.get_mp()) ;
154 
155  double erreur = 0 ;
156  for (int comp=0 ; comp<3 ; comp++)
157  {
158  val_hor.annule_hard() ; // Pour initialiser les tableaux
159  for (int k=0 ; k<np ; k++)
160  for (int j=0 ; j<nt ; j++)
161  for (int i=0 ; i<nr ; i++)
162  val_hor.set(1, k, j, i) = derive_r(comp) (1, k, j, 0) ;
163 
164  double r_0 = shift_auto.get_mp()->val_r (1, -1, 0, 0) ;
165  double r_1 = shift_auto.get_mp()->val_r (1, 1, 0, 0) ;
166 
167  fonction_radiale = pow(r_1-shift_auto.get_mp()->r, 3.)*
168  (shift_auto.get_mp()->r-r_0)/pow(r_1-r_0, 3.) ;
169  fonction_radiale.annule(0) ;
170  fonction_radiale.annule(2, nz-1) ;
171 
172  enleve = fonction_radiale * val_hor ;
173  enleve.va.base = shift_auto(comp).va.base ;
174 
175  if (norme(enleve)(1) != 0)
176  shift_auto.set(comp) = shift_auto(comp) - enleve ;
177  if (norme(shift_auto(comp))(1) > 1e-5) {
178  Tbl diff (diffrelmax (shift_auto(comp), shift_old(comp))) ;
179  if (erreur < diff(1))
180  erreur = diff(1) ;
181  }
182  }
183  return erreur ;
184 }
185 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:539
Lorene::norme
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:481
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Cmp::annule
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348