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
23char 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
73namespace Lorene {
74void 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}
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:539
Lorene prototypes.
Definition app_hor.h:64