source: src/discretization/discretization_poisson_fv.cpp@ 14d38c

Last change on this file since 14d38c was 14d38c, checked in by Julian Iseringhausen <isering@…>, 13 years ago

Renamed folder samples to discretization and worked on inner boundaries for open boundary conditions.
  • Property mode set to 100644
File size: 15.5 KB
Line 
1/*
2 * vmg - a versatile multigrid solver
3 * Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
4 *
5 * vmg is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * vmg is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19/**
20 * @file discretization_poisson_fv.cpp
21 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
22 * @date Mon Apr 18 13:03:47 2011
23 *
24 * @brief Finite volume discretization for the Poisson
25 * equation. Absolutely equivalent to the finite
26 * difference discretization unless you use
27 * hierarchically coarsened grids.
28 *
29 */
30
31#ifdef HAVE_CONFIG_H
32#include <config.h>
33#endif
34
35#include "comm/comm.hpp"
36#include "discretization/discretization_poisson_fv.hpp"
37#include "grid/grid_index_translations.hpp"
38
39using namespace VMG;
40
41void DiscretizationPoissonFV::InitDiscretizationPoissonFV()
42{
43 switch (order)
44 {
45 case 2:
46 stencil.SetDiag(6.0);
47 stencil.push_back(-1, 0, 0, -1.0);
48 stencil.push_back( 1, 0, 0, -1.0);
49 stencil.push_back( 0, -1, 0, -1.0);
50 stencil.push_back( 0, 1, 0, -1.0);
51 stencil.push_back( 0, 0, -1, -1.0);
52 stencil.push_back( 0, 0, 1, -1.0);
53 break;
54 case 4:
55 stencil.SetDiag(24.0/6.0);
56 stencil.push_back(-1, 0, 0, -2.0/6.0);
57 stencil.push_back( 1, 0, 0, -2.0/6.0);
58 stencil.push_back( 0, -1, 0, -2.0/6.0);
59 stencil.push_back( 0, 1, 0, -2.0/6.0);
60 stencil.push_back( 0, 0, -1, -2.0/6.0);
61 stencil.push_back( 0, 0, 1, -2.0/6.0);
62 stencil.push_back(-1, -1, 0, -1.0/6.0);
63 stencil.push_back(-1, 1, 0, -1.0/6.0);
64 stencil.push_back( 1, -1, 0, -1.0/6.0);
65 stencil.push_back( 1, 1, 0, -1.0/6.0);
66 stencil.push_back(-1, 0, -1, -1.0/6.0);
67 stencil.push_back(-1, 0, 1, -1.0/6.0);
68 stencil.push_back( 1, 0, -1, -1.0/6.0);
69 stencil.push_back( 1, 0, 1, -1.0/6.0);
70 stencil.push_back( 0, -1, -1, -1.0/6.0);
71 stencil.push_back( 0, -1, 1, -1.0/6.0);
72 stencil.push_back( 0, 1, -1, -1.0/6.0);
73 stencil.push_back( 0, 1, 1, -1.0/6.0);
74 break;
75 default:
76 assert(0 != "vmg choose discretization order 2 or 4");
77 break;
78 }
79}
80
81void DiscretizationPoissonFV::ModifyRightHandSide()
82{
83 if (order == 4) {
84
85 Grid& rhs = MG::GetRhsMaxLevel();
86
87 Stencil stencil(6.0/12.0);
88 stencil.push_back(-1, 0, 0, 1.0/12.0);
89 stencil.push_back( 1, 0, 0, 1.0/12.0);
90 stencil.push_back( 0, -1, 0, 1.0/12.0);
91 stencil.push_back( 0, 1, 0, 1.0/12.0);
92 stencil.push_back( 0, 0, -1, 1.0/12.0);
93 stencil.push_back( 0, 0, 1, 1.0/12.0);
94
95 stencil.Apply(rhs);
96
97 }
98}
99
100void DiscretizationPoissonFV::SetInnerBoundaryCompute(Grid& sol_f, Grid& rhs_f, Grid& sol_c) const
101{
102 Index i_c, i_f;
103
104 const Boundary& bc = MG::GetComm()->BoundaryConditions();
105 const Index off((GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[0] % 2 == 0 ? 0 : 1),
106 (GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[1] % 2 == 0 ? 0 : 1),
107 (GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[2] % 2 == 0 ? 0 : 1));
108
109 const Index begin_f = sol_f.Local().Begin() + off;
110 const Index end_f = sol_f.Local().End();
111
112 const Index begin_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, GridIndexTranslations::LocalToGlobalFinest(sol_f, begin_f));
113
114 const Index b1_finest = GridIndexTranslations::LocalToGlobalFinest(sol_f, sol_f.Local().Begin());
115 const Index b2_finest = GridIndexTranslations::LocalToGlobalFinest(sol_f, sol_f.Local().End()-1);
116
117 const Index b1_f = sol_f.Local().BoundaryBegin1();
118 const Index b2_f = sol_f.Local().BoundaryBegin2();
119
120 const Index b1_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, b1_finest) - 1;
121 const Index b2_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, b2_finest) + 1;
122
123 const vmg_float c_1_3 = 1.0 / 3.0;
124 const vmg_float c_2_3 = 2.0 / 3.0;
125 const vmg_float c_4_3 = 4.0 / 3.0;
126
127
128 //
129 // X-direction
130 //
131 if (bc.X() == Open) {
132
133 if (sol_f.Local().BoundarySize1().X() > 0) {
134
135 for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
136 for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
137 rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (sol_c(b1_c.X()-1,i_c.Y(),i_c.Z()) - sol_f(b1_f.X()+1,i_f.Y(),i_f.Z()));
138
139 for (i_f.Y()=begin_f.Y()-1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
140 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
141 rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()) + rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()));
142
143 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
144 for (i_f.Z()=begin_f.Z()-1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
145 rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b1_f.X(),i_f.Y(),i_f.Z()-1) + rhs_f(b1_f.X(),i_f.Y(),i_f.Z()+1));
146
147 for (i_f.Y()=begin_f.Y()-1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
148 for (i_f.Z()=begin_f.Z()-1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
149 rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.25 * (rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()-1) +
150 rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()-1) +
151 rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()+1) +
152 rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()+1));
153
154 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
155 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
156 rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = sol_f(b1_f.X(),i_f.Y(),i_f.Z()) = c_2_3 * rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) +
157 c_4_3 * sol_f(b1_f.X()+1,i_f.Y(),i_f.Z()) -
158 c_1_3 * sol_f(b1_f.X()+2,i_f.Y(),i_f.Z());
159
160 }
161
162 if (sol_f.Local().BoundarySize2().X() > 0) {
163
164 for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
165 for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
166 rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (sol_c(b2_c.X()+1,i_c.Y(),i_c.Z()) - sol_f(b2_f.X()-1,i_f.Y(),i_f.Z()));
167
168 for (i_f.Y()=begin_f.Y()-1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
169 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
170 rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()) + rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()));
171
172 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
173 for (i_f.Z()=begin_f.Z()-1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
174 rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b2_f.X(),i_f.Y(),i_f.Z()-1) + rhs_f(b2_f.X(),i_f.Y(),i_f.Z()+1));
175
176 for (i_f.Y()=begin_f.Y()-1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
177 for (i_f.Z()=begin_f.Z()-1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
178 rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.25 * (rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()-1) +
179 rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()-1) +
180 rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()+1) +
181 rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()+1));
182
183 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
184 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
185 rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = sol_f(b2_f.X(),i_f.Y(),i_f.Z()) = c_2_3 * rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) +
186 c_4_3 * sol_f(b2_f.X()-1,i_f.Y(),i_f.Z()) -
187 c_1_3 * sol_f(b2_f.X()-2,i_f.Y(),i_f.Z());
188
189 }
190 }
191
192 //
193 // Y-direction
194 //
195 if (bc.Y() == Open) {
196
197 if (sol_f.Local().BoundarySize1().Y() > 0) {
198
199 for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
200 for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
201 rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (sol_c(i_c.X(),b1_c.Y()-1,i_c.Z()) - sol_f(i_f.X(),b1_f.Y()+1,i_f.Z()));
202
203 for (i_f.X()=begin_f.X()-1; i_f.X()<end_f.X(); i_f.X()+=2)
204 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
205 rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()) + rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()));
206
207 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
208 for (i_f.Z()=begin_f.Z()-1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
209 rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X(),b1_f.Y(),i_f.Z()-1) + rhs_f(i_f.X(),b1_f.Y(),i_f.Z()+1));
210
211 for (i_f.X()=begin_f.X()-1; i_f.X()<end_f.X(); i_f.X()+=2)
212 for (i_f.Z()=begin_f.Z()-1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
213 rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()-1) +
214 rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()-1) +
215 rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()+1) +
216 rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()+1));
217
218 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
219 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
220 rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = sol_f(i_f.X(),b1_f.Y(),i_f.Z()) = c_2_3 * rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) +
221 c_4_3 * sol_f(i_f.X(),b1_f.Y()+1,i_f.Z()) -
222 c_1_3 * sol_f(i_f.X(),b1_f.Y()+2,i_f.Z());
223
224 }
225
226 if (sol_f.Local().BoundarySize2().Y() > 0) {
227
228 for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
229 for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
230 rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (sol_c(i_c.X(),b2_c.Y()+1,i_c.Z()) - sol_f(i_f.X(),b2_f.Y()-1,i_f.Z()));
231
232 for (i_f.X()=begin_f.X()-1; i_f.X()<end_f.X(); i_f.X()+=2)
233 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
234 rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()) + rhs_f(i_f.X()+1,b2_f.Z(),i_f.Z()));
235
236 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
237 for (i_f.Z()=begin_f.Z()-1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
238 rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X(),b2_f.Y(),i_f.Z()-1) + rhs_f(i_f.X(),b2_f.Y(),i_f.Z()+1));
239
240 for (i_f.X()=begin_f.X()-1; i_f.X()<end_f.X(); i_f.X()+=2)
241 for (i_f.Z()=begin_f.Z()-1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
242 rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()-1) +
243 rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()-1) +
244 rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()+1) +
245 rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()+1));
246
247 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
248 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
249 rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = sol_f(i_f.X(),b2_f.Y(),i_f.Z()) = c_2_3 * rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) +
250 c_4_3 * sol_f(i_f.X(),b2_f.Y()-1,i_f.Z()) -
251 c_1_3 * sol_f(i_f.X(),b2_f.Y()-2,i_f.Z());
252
253 }
254
255 }
256
257 //
258 // Z-direction
259 //
260 if (bc.Z() == Open) {
261
262 if (sol_f.Local().BoundarySize1().Z() > 0) {
263
264 for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
265 for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
266 rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (sol_c(i_c.X(),i_c.Y(),b1_c.Z()-1) - sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1));
267
268 for (i_f.X()=begin_f.X()-1; i_f.X()<end_f.X(); i_f.X()+=2)
269 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
270 rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,i_f.Y(),b1_f.Z()) + rhs_f(i_f.X()+1,i_f.Y(),b1_f.Z()));
271
272 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
273 for (i_f.Y()=begin_f.Y()-1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
274 rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (rhs_f(i_f.X(),i_f.Y()-1,b1_f.Z()) + rhs_f(i_f.X(),i_f.Y()+1,b1_f.Z()));
275
276 for (i_f.X()=begin_f.X()-1; i_f.X()<end_f.X(); i_f.X()+=2)
277 for (i_f.Y()=begin_f.Y()-1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
278 rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,i_f.Y()-1,b1_f.Z()) +
279 rhs_f(i_f.X()+1,i_f.Y()-1,b1_f.Z()) +
280 rhs_f(i_f.X()-1,i_f.Y()+1,b1_f.Z()) +
281 rhs_f(i_f.X()+1,i_f.Y()+1,b1_f.Z()));
282
283 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
284 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
285 rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = sol_f(i_f.X(),i_f.Y(),b1_f.Z()) = c_2_3 * rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) +
286 c_4_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1) -
287 c_1_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+2);
288
289 }
290
291 if (sol_f.Local().BoundarySize2().Z() > 0) {
292
293 for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
294 for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
295 rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (sol_c(i_c.X(),i_c.Y(),b2_c.Z()+1) - sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1));
296
297 for (i_f.X()=begin_f.X()-1; i_f.X()<end_f.X(); i_f.X()+=2)
298 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
299 rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,i_f.Y(),b2_f.Z()) + rhs_f(i_f.X()+1,i_f.Y(),b2_f.Z()));
300
301 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
302 for (i_f.Y()=begin_f.Y()-1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
303 rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (rhs_f(i_f.X(),i_f.Y()-1,b2_f.Z()) + rhs_f(i_f.X(),i_f.Y()+1,b2_f.Z()));
304
305 for (i_f.X()=begin_f.X()-1; i_f.X()<end_f.X(); i_f.X()+=2)
306 for (i_f.Y()=begin_f.Y()-1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
307 rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,i_f.Y()-1,b2_f.Z()) +
308 rhs_f(i_f.X()+1,i_f.Y()-1,b2_f.Z()) +
309 rhs_f(i_f.X()-1,i_f.Y()+1,b2_f.Z()) +
310 rhs_f(i_f.X()+1,i_f.Y()+1,b2_f.Z()));
311
312 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
313 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
314 rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = sol_f(i_f.X(),i_f.Y(),b2_f.Z()) = c_2_3 * rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) +
315 c_4_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1) -
316 c_1_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-2);
317
318 }
319
320 }
321
322#ifdef DEBUG_MATRIX_CHECKS
323 rhs_f.IsConsistent();
324 sol_f.IsConsistent();
325#endif
326}
Note: See TracBrowser for help on using the repository browser.