source: ThirdParty/vmg/src/discretization/discretization_poisson_fv.cpp@ 33a694

Action_Thermostats Add_AtomRandomPerturbation Add_RotateAroundBondAction Add_SelectAtomByNameAction Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_StructOpt_integration_tests AutomationFragmentation_failures Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator Docu_Python_wait EmpiricalPotential_contain_HomologyGraph_documentation Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph Fix_ChronosMutex Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion GeometryObjects Gui_displays_atomic_force_velocity IndependentFragmentGrids_IntegrationTest JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks RotateToPrincipalAxisSystem_UndoRedo StoppableMakroAction TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps Ubuntu_1604_changes
Last change on this file since 33a694 was 7faa5c, checked in by Frederik Heber <heber@…>, 8 years ago

Merge commit 'de061d9d851257a04e924d4472df4523d33bb08b' as 'ThirdParty/vmg'
  • Property mode set to 100644
File size: 15.6 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 <libvmg_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 Comm& comm = *MG::GetComm();
105
106 const Boundary& bc = comm.BoundaryConditions();
107 const Index off((GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[0] % 2 == 0 ? 0 : 1),
108 (GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[1] % 2 == 0 ? 0 : 1),
109 (GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[2] % 2 == 0 ? 0 : 1));
110
111 const Index begin_f = sol_f.Local().Begin() - off;
112 const Index end_f = sol_f.Local().End();
113
114 const Index begin_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, GridIndexTranslations::LocalToGlobalFinest(sol_f, begin_f));
115
116 const Index b1_f = sol_f.Local().BoundaryBegin1();
117 const Index b2_f = sol_f.Local().BoundaryBegin2();
118
119 const Index b1_finest = GridIndexTranslations::LocalToGlobalFinest(sol_f, sol_f.Local().Begin() + off);
120 const Index b2_finest = GridIndexTranslations::LocalToGlobalFinest(sol_f, sol_f.Local().End()-1);
121
122 const Index b1_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, b1_finest) - 1;
123 const Index b2_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, b2_finest) + 1;
124
125 const vmg_float c_1_3 = 1.0 / 3.0;
126 const vmg_float c_2_3 = 2.0 / 3.0;
127 const vmg_float c_4_3 = 4.0 / 3.0;
128
129 comm.CommToGhosts(sol_f);
130 comm.CommToGhosts(sol_c);
131
132 //
133 // X-direction
134 //
135 if (bc.X() == Open) {
136
137 if (sol_f.Local().BoundarySize1().X() > 0) {
138
139 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())
140 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())
141 rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (sol_c(b1_c.X(),i_c.Y(),i_c.Z()) - sol_f(b1_f.X()+1,i_f.Y(),i_f.Z()));
142
143 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
144 for (i_f.Z()=begin_f.Z(); 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()-1,i_f.Z()) + rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()));
146
147 for (i_f.Y()=begin_f.Y(); 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.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));
150
151 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
152 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
153 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) +
154 rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()-1) +
155 rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()+1) +
156 rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()+1));
157
158 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
159 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
160 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()) +
161 c_4_3 * sol_f(b1_f.X()+1,i_f.Y(),i_f.Z()) -
162 c_1_3 * sol_f(b1_f.X()+2,i_f.Y(),i_f.Z());
163
164 }
165
166 if (sol_f.Local().BoundarySize2().X() > 0) {
167
168 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())
169 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())
170 rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (sol_c(b2_c.X(),i_c.Y(),i_c.Z()) - sol_f(b2_f.X()-1,i_f.Y(),i_f.Z()));
171
172 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
173 for (i_f.Z()=begin_f.Z(); 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()-1,i_f.Z()) + rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()));
175
176 for (i_f.Y()=begin_f.Y(); 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.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));
179
180 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
181 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
182 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) +
183 rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()-1) +
184 rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()+1) +
185 rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()+1));
186
187 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
188 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
189 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()) +
190 c_4_3 * sol_f(b2_f.X()-1,i_f.Y(),i_f.Z()) -
191 c_1_3 * sol_f(b2_f.X()-2,i_f.Y(),i_f.Z());
192
193 }
194 }
195
196 //
197 // Y-direction
198 //
199 if (bc.Y() == Open) {
200
201 if (sol_f.Local().BoundarySize1().Y() > 0) {
202
203 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())
204 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())
205 rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (sol_c(i_c.X(),b1_c.Y(),i_c.Z()) - sol_f(i_f.X(),b1_f.Y()+1,i_f.Z()));
206
207 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
208 for (i_f.Z()=begin_f.Z(); 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()-1,b1_f.Y(),i_f.Z()) + rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()));
210
211 for (i_f.X()=begin_f.X(); 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.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));
214
215 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
216 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
217 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) +
218 rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()-1) +
219 rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()+1) +
220 rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()+1));
221
222 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
223 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
224 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()) +
225 c_4_3 * sol_f(i_f.X(),b1_f.Y()+1,i_f.Z()) -
226 c_1_3 * sol_f(i_f.X(),b1_f.Y()+2,i_f.Z());
227
228 }
229
230 if (sol_f.Local().BoundarySize2().Y() > 0) {
231
232 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())
233 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())
234 rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (sol_c(i_c.X(),b2_c.Y(),i_c.Z()) - sol_f(i_f.X(),b2_f.Y()-1,i_f.Z()));
235
236 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
237 for (i_f.Z()=begin_f.Z(); 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()-1,b2_f.Y(),i_f.Z()) + rhs_f(i_f.X()+1,b2_f.Z(),i_f.Z()));
239
240 for (i_f.X()=begin_f.X(); 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.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));
243
244 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
245 for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
246 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) +
247 rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()-1) +
248 rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()+1) +
249 rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()+1));
250
251 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
252 for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
253 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()) +
254 c_4_3 * sol_f(i_f.X(),b2_f.Y()-1,i_f.Z()) -
255 c_1_3 * sol_f(i_f.X(),b2_f.Y()-2,i_f.Z());
256
257 }
258
259 }
260
261 //
262 // Z-direction
263 //
264 if (bc.Z() == Open) {
265
266 if (sol_f.Local().BoundarySize1().Z() > 0) {
267
268 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())
269 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())
270 rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (sol_c(i_c.X(),i_c.Y(),b1_c.Z()) - sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1));
271
272 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
273 for (i_f.Y()=begin_f.Y(); 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()-1,i_f.Y(),b1_f.Z()) + rhs_f(i_f.X()+1,i_f.Y(),b1_f.Z()));
275
276 for (i_f.X()=begin_f.X(); 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.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()));
279
280 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
281 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
282 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()) +
283 rhs_f(i_f.X()+1,i_f.Y()-1,b1_f.Z()) +
284 rhs_f(i_f.X()-1,i_f.Y()+1,b1_f.Z()) +
285 rhs_f(i_f.X()+1,i_f.Y()+1,b1_f.Z()));
286
287 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
288 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
289 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()) +
290 c_4_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1) -
291 c_1_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+2);
292
293 }
294
295 if (sol_f.Local().BoundarySize2().Z() > 0) {
296
297 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())
298 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())
299 rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (sol_c(i_c.X(),i_c.Y(),b2_c.Z()) - sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1));
300
301 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
302 for (i_f.Y()=begin_f.Y(); 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()-1,i_f.Y(),b2_f.Z()) + rhs_f(i_f.X()+1,i_f.Y(),b2_f.Z()));
304
305 for (i_f.X()=begin_f.X(); 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.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()));
308
309 for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
310 for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
311 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()) +
312 rhs_f(i_f.X()+1,i_f.Y()-1,b2_f.Z()) +
313 rhs_f(i_f.X()-1,i_f.Y()+1,b2_f.Z()) +
314 rhs_f(i_f.X()+1,i_f.Y()+1,b2_f.Z()));
315
316 for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
317 for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
318 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()) +
319 c_4_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1) -
320 c_1_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-2);
321
322 }
323
324 }
325
326#ifdef DEBUG_MATRIX_CHECKS
327 rhs_f.IsConsistent();
328 sol_f.IsConsistent();
329#endif
330}
Note: See TracBrowser for help on using the repository browser.