17 #ifndef __deal2__multigrid_templates_h
18 #define __deal2__multigrid_templates_h
19 #include <deal.II/multigrid/multigrid.h>
21 #include <deal.II/base/logstream.h>
25 DEAL_II_NAMESPACE_OPEN
28 template <
class VECTOR>
30 const unsigned int maxlevel,
41 defect(minlevel,maxlevel),
42 solution(minlevel,maxlevel),
44 matrix(&matrix, typeid(*this).name()),
45 coarse(&coarse, typeid(*this).name()),
46 transfer(&transfer, typeid(*this).name()),
47 pre_smooth(&pre_smooth, typeid(*this).name()),
48 post_smooth(&post_smooth, typeid(*this).name()),
49 edge_out(0, typeid(*this).name()),
50 edge_in(0, typeid(*this).name()),
51 edge_down(0, typeid(*this).name()),
52 edge_up(0, typeid(*this).name()),
58 template <
class VECTOR>
61 const unsigned int max_level)
65 defect.resize(minlevel, maxlevel);
69 template <
class VECTOR>
79 template <
class VECTOR>
91 template <
class VECTOR>
99 template <
class VECTOR>
107 template <
class VECTOR>
117 template <
class VECTOR>
127 template <
class VECTOR>
132 deallog <<
"V-cycle entering level " << level << std::endl;
134 deallog <<
"V-cycle Defect norm " << defect[level].l2_norm()
137 if (level == minlevel)
140 deallog <<
"Coarse level " << level << std::endl;
141 (*coarse)(level, solution[level], defect[level]);
145 deallog <<
"Smoothing on level " << level << std::endl;
150 pre_smooth->smooth(level, solution[level], defect[level]);
154 deallog <<
"Solution norm " << solution[level].l2_norm()
158 deallog <<
"Residual on level " << level << std::endl;
160 matrix->vmult(level, t[level], solution[level]);
163 deallog <<
"Residual norm " << t[level].l2_norm()
173 for (
unsigned int l = level; l>minlevel; --l)
176 if (l==level && edge_out != 0)
178 edge_out->vmult_add(level, t[level], solution[level]);
180 deallog <<
"Norm t[" << level <<
"] " << t[level].l2_norm() << std::endl;
183 if (l==level && edge_down != 0)
184 edge_down->vmult(level, t[level-1], solution[level]);
186 transfer->restrict_and_add (l, t[l-1], t[l]);
188 deallog <<
"restrict t[" << l-1 <<
"] " << t[l-1].l2_norm() << std::endl;
189 defect[l-1] -= t[l-1];
191 deallog <<
"defect d[" << l-1 <<
"] " << defect[l-1].l2_norm() << std::endl;
195 solution[level-1] = 0.;
196 level_v_step(level-1);
205 transfer->prolongate(level, t[level], solution[level-1]);
207 deallog <<
"Prolongate norm " << t[level].l2_norm() << std::endl;
209 solution[level] += t[level];
213 edge_in->Tvmult(level, t[level], solution[level]);
214 defect[level] -= t[level];
219 edge_up->Tvmult(level, t[level], solution[level-1]);
220 defect[level] -= t[level];
224 deallog <<
"V-cycle Defect norm " << defect[level].l2_norm()
228 deallog <<
"Smoothing on level " << level << std::endl;
233 post_smooth->smooth(level, solution[level], defect[level]);
238 deallog <<
"Solution norm " << solution[level].l2_norm()
242 deallog <<
"V-cycle leaving level " << level << std::endl;
247 template <
class VECTOR>
267 deallog << cychar <<
"-cycle entering level " << level << std::endl;
274 defect2[level-1] = 0.;
275 transfer->restrict_and_add (level, defect2[level-1], defect2[level]);
283 t[level].equ(-1.,defect2[level],1.,defect[level]);
286 deallog << cychar <<
"-cycle defect norm " << t[level].l2_norm()
289 if (level == minlevel)
292 deallog << cychar <<
"-cycle coarse level " << level << std::endl;
294 (*coarse)(level, solution[level], t[level]);
298 deallog << cychar <<
"-cycle smoothing level " << level << std::endl;
301 pre_smooth->smooth(level, solution[level], t[level]);
304 deallog << cychar <<
"-cycle solution norm "
305 << solution[level].l2_norm() << std::endl;
308 deallog << cychar <<
"-cycle residual level " << level << std::endl;
310 matrix->vmult(level, t[level], solution[level]);
317 edge_out->vmult_add(level, t[level], solution[level]);
320 edge_down->vmult_add(level, defect2[level-1], solution[level]);
322 transfer->restrict_and_add (level, defect2[level-1], t[level]);
324 solution[level-1] = 0.;
329 level_step(level-1, cycle);
333 if (level>minlevel+1)
336 if (cycle == w_cycle)
337 level_step(level-1, cycle);
340 else if (cycle == f_cycle)
341 level_step(level-1, v_cycle);
350 transfer->prolongate(level, t[level], solution[level-1]);
352 solution[level] += t[level];
356 edge_in->Tvmult(level, t[level], solution[level]);
359 edge_up->Tvmult(level, t[level], solution[level-1]);
361 t[level].sadd(-1.,-1.,defect2[level],1.,defect[level]);
364 deallog << cychar <<
"-cycle Defect norm " << t[level].l2_norm()
368 deallog << cychar <<
"-cycle smoothing level " << level << std::endl;
370 post_smooth->smooth(level, solution[level], t[level]);
373 deallog << cychar <<
"-cycle leaving level " << level << std::endl;
377 template <
class VECTOR>
384 solution.resize(minlevel, maxlevel);
385 t.resize(minlevel, maxlevel);
386 if (cycle_type != v_cycle)
387 defect2.resize(minlevel, maxlevel);
389 for (
unsigned int level=minlevel; level<=maxlevel; ++level)
391 solution[level].reinit(defect[level]);
392 t[level].reinit(defect[level]);
393 if (cycle_type != v_cycle)
394 defect2[level].reinit(defect[level]);
397 if (cycle_type == v_cycle)
398 level_v_step (maxlevel);
400 level_step (maxlevel, cycle_type);
404 template <
class VECTOR>
411 solution.resize(minlevel, maxlevel);
412 t.resize(minlevel, maxlevel);
414 for (
unsigned int level=minlevel; level<=maxlevel; ++level)
416 solution[level].reinit(defect[level]);
417 t[level].reinit(defect[level]);
419 level_v_step (maxlevel);
423 template <
class VECTOR>
428 mg.
defect[maxlevel] = src;
429 for (
unsigned int level=maxlevel; level>minlevel; --level)
432 mg.
transfer->restrict_and_add (level,
442 template <
class VECTOR>
447 mg.
defect[maxlevel] = src;
453 template <
class VECTOR>
461 template <
class VECTOR>
468 DEAL_II_NAMESPACE_CLOSE
Multigrid(const DoFHandler< dim > &mg_dof_handler, const MGMatrixBase< VECTOR > &matrix, const MGCoarseGridBase< VECTOR > &coarse, const MGTransferBase< VECTOR > &transfer, const MGSmootherBase< VECTOR > &pre_smooth, const MGSmootherBase< VECTOR > &post_smooth, Cycle cycle=v_cycle)
void Tvmult_add(VECTOR &dst, const VECTOR &src) const DEAL_II_DEPRECATED
void set_edge_matrices(const MGMatrixBase< VECTOR > &edge_out, const MGMatrixBase< VECTOR > &edge_in)
MGLevelObject< VECTOR > solution
void Tvmult(VECTOR &dst, const VECTOR &src) const DEAL_II_DEPRECATED
void level_step(const unsigned int level, Cycle cycle)
void level_v_step(const unsigned int level)
void vmult_add(VECTOR &dst, const VECTOR &src) const DEAL_II_DEPRECATED
SmartPointer< const MGTransferBase< VECTOR >, Multigrid< VECTOR > > transfer
#define Assert(cond, exc)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void vmult(VECTOR &dst, const VECTOR &src) const DEAL_II_DEPRECATED
void set_maxlevel(const unsigned int)
void set_minlevel(const unsigned int level, bool relative=false)
void reinit(const unsigned int minlevel, const unsigned int maxlevel)
void set_edge_flux_matrices(const MGMatrixBase< VECTOR > &edge_down, const MGMatrixBase< VECTOR > &edge_up)
::ExceptionBase & ExcNotImplemented()
void set_debug(const unsigned int)
MGLevelObject< VECTOR > defect