LORENE
separation.C
1 /*
2  * Copyright (c) 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 separation_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/separation.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $" ;
24 
25 /*
26  * $Id: separation.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $
27  * $Log: separation.C,v $
28  * Revision 1.4 2014/10/13 08:52:41 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.3 2014/10/06 15:12:59 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.2 2003/10/03 15:58:44 j_novak
35  * Cleaning of some headers
36  *
37  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
38  * LORENE
39  *
40  * Revision 2.6 2001/04/02 12:16:20 phil
41  * *** empty log message ***
42  *
43  * Revision 2.5 2001/03/30 13:48:17 phil
44  * on appelle raccord externe
45  *
46  * Revision 2.4 2001/03/22 10:40:30 phil
47  * modification prototypage
48  *
49  * Revision 2.3 2001/03/02 10:19:05 phil
50  * modification parametrage pour affichage
51  *
52  * Revision 2.2 2001/02/28 13:39:34 phil
53  * modif cas etat_zero
54  *
55  * Revision 2.1 2001/02/28 13:23:00 phil
56  * modif etat initial
57  *
58  * Revision 2.0 2001/02/28 11:24:34 phil
59  * *** empty log message ***
60  *
61  *
62  * $Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/separation.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $
63  *
64  */
65 
66 //standard
67 #include <cstdlib>
68 
69 // Lorene
70 #include "cmp.h"
71 #include "proto.h"
72 
73 namespace Lorene {
74 void separation (const Cmp& c1, const Cmp& c2, Cmp& res1, Cmp& res2, int decrois,
75  int puiss, int lmax, double precision, const double relax, const int itemax, const int flag) {
76 
77  assert (c1.get_etat() != ETATNONDEF) ;
78  assert (c2.get_etat() != ETATNONDEF) ;
79 
80  if ((c1.get_etat() == ETATZERO) && (c2.get_etat() == ETATZERO)) {
81  res1.set_etat_zero() ;
82  res2.set_etat_zero() ;
83  return ;
84  }
85  else {
86 
87  res1 = c1 ;
88  if (res1.get_etat() == ETATZERO) {
89  res1.annule_hard() ;
90  res1.std_base_scal() ;
91  }
92 
93  res1.raccord_externe (decrois, puiss, lmax) ;
94  for (int i=0 ; i<decrois ; i++)
95  res1.dec_dzpuis() ;
96 
97  res2 = c2 ;
98  if (res2.get_etat() == ETATZERO) {
99  res2.annule_hard() ;
100  res2.std_base_scal() ;
101  }
102  res2.raccord_externe (decrois, puiss, lmax) ;
103  for (int i=0 ; i<decrois ; i++)
104  res2.dec_dzpuis() ;
105 
106  int indic = 1 ;
107  int conte = 0 ;
108  // On commence la boucle pour separer :
109  while (indic == 1) {
110  Cmp old_un (res1) ;
111  Cmp old_deux (res2) ;
112 
113  // On fait les modifications :
114  Mtbl xa_mtbl_un (c1.get_mp()->xa) ;
115  Mtbl ya_mtbl_un (c1.get_mp()->ya) ;
116  Mtbl za_mtbl_un (c1.get_mp()->za) ;
117 
118  Mtbl xa_mtbl_deux (c2.get_mp()->xa) ;
119  Mtbl ya_mtbl_deux (c2.get_mp()->ya) ;
120  Mtbl za_mtbl_deux (c2.get_mp()->za) ;
121 
122  double xabs, yabs, zabs, air, theta, phi ;
123  int np, nt, nr ;
124 
125  // On modifie le Cmp 1
126  int nz_un = c1.get_mp()->get_mg()->get_nzone() ;
127  for (int l=1 ; l<nz_un-1 ; l++) {
128 
129  np = c1.get_mp()->get_mg()->get_np(l) ;
130  nt = c1.get_mp()->get_mg()->get_nt(l) ;
131  nr = c1.get_mp()->get_mg()->get_nr(l) ;
132 
133  for (int k=0 ; k<np ; k++)
134  for (int j=0 ; j<nt ; j++)
135  for (int i=0 ; i<nr ; i++) {
136  xabs = xa_mtbl_un (l, k, j, i) ;
137  yabs = ya_mtbl_un (l, k, j, i) ;
138  zabs = za_mtbl_un (l, k, j, i) ;
139 
140  c2.get_mp()->convert_absolute(xabs, yabs, zabs, air, theta, phi) ;
141  res1.set(l, k, j, i) =
142  (1-relax)*res1.set(l, k, j, i) +
143  relax*(c1(l, k, j, i) - old_deux.val_point(air, theta, phi)) ;
144  }
145  }
146 
147  // On modifie le trou 2
148  int nz_deux = c2.get_mp()->get_mg()->get_nzone() ;
149  for (int l=1 ; l<nz_deux-1 ; l++) {
150 
151  np = c2.get_mp()->get_mg()->get_np(l) ;
152  nt = c2.get_mp()->get_mg()->get_nt(l) ;
153  nr = c2.get_mp()->get_mg()->get_nr(l) ;
154 
155  for (int k=0 ; k<np ; k++)
156  for (int j=0 ; j<nt ; j++)
157  for (int i=0 ; i<nr ; i++) {
158 
159  xabs = xa_mtbl_deux (l, k, j, i) ;
160  yabs = ya_mtbl_deux (l, k, j, i) ;
161  zabs = za_mtbl_deux (l, k, j, i) ;
162 
163  c1.get_mp()->convert_absolute(xabs, yabs, zabs, air, theta, phi) ;
164  res2.set(l, k, j, i) =
165  (1-relax)*res2.set(l, k, j, i) +
166  relax*(c2(l, k, j, i) - old_un.val_point(air, theta, phi)) ;
167  }
168  }
169 
170  // les coefficients ne sont plus a jour :
171  res1.va.set_etat_c_qcq() ;
172  res2.va.set_etat_c_qcq() ;
173  // On raccord dans la zec :
174  res1.raccord_externe (decrois, puiss, lmax) ;
175  for (int i=0 ; i<decrois ; i++)
176  res1.dec_dzpuis() ;
177 
178  res1.va.coef_i() ;
179  res2.raccord_externe (decrois, puiss, lmax) ;
180  for (int i=0 ; i<decrois ; i++)
181  res2.dec_dzpuis() ;
182  res2.va.coef_i() ;
183 
184  // On regarde si on a converge :
185 
186  double erreur = 0 ;
187 
188  Tbl diff_un (diffrelmax(res1, old_un)) ;
189  for (int i=1 ; i<nz_un-1 ; i++)
190  if (diff_un(i)>erreur)
191  erreur = diff_un(i) ;
192 
193  Tbl diff_deux (diffrelmax(res2, old_deux)) ;
194  for (int i=1 ; i<nz_deux-1 ; i++)
195  if (diff_deux(i)>erreur)
196  erreur = diff_deux(i) ;
197 
198  if (flag == 1)
199  cout << "Pas " << conte << " : erreur = " << erreur << endl ;
200  if (erreur<=precision)
201  indic = -1 ;
202 
203  conte ++ ;
204  if (conte > itemax)
205  indic = -1 ;
206  }
207  }
208 }
209 }
Lorene prototypes.
Definition: app_hor.h:64
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:539