1 /* ---------------------------------------------------------------------
2 *
3 * Copyright (C) 2007 - 2020 by the deal.II authors
4 *
5 * This file is part of the deal.II library.
6 *
7 * The deal.II library is free software; you can use it, redistribute
8 * it, and/or modify it under the terms of the GNU Lesser General
9 * Public License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 * The full text of the license can be found in the file LICENSE.md at
12 * the top level directory of deal.II.
13 *
14 * ---------------------------------------------------------------------
15
16 *
17 * Author: Tobias Leicht, 2007
18 */
19
20
21 // The deal.II include files have already been covered in previous examples
22 // and will thus not be further commented on.
23 #include <deal.II/base/function.h>
24 #include <deal.II/base/quadrature_lib.h>
25 #include <deal.II/base/timer.h>
26 #include <deal.II/lac/precondition_block.h>
27 #include <deal.II/lac/solver_richardson.h>
28 #include <deal.II/lac/sparse_matrix.h>
29 #include <deal.II/lac/vector.h>
30 #include <deal.II/grid/tria.h>
31 #include <deal.II/grid/grid_generator.h>
32 #include <deal.II/grid/grid_out.h>
33 #include <deal.II/grid/grid_refinement.h>
34 #include <deal.II/grid/tria_accessor.h>
35 #include <deal.II/grid/tria_iterator.h>
36 #include <deal.II/dofs/dof_handler.h>
37 #include <deal.II/dofs/dof_accessor.h>
38 #include <deal.II/dofs/dof_tools.h>
39 #include <deal.II/fe/fe_values.h>
40 #include <deal.II/fe/mapping_q1.h>
41 #include <deal.II/fe/fe_dgq.h>
42 #include <deal.II/numerics/data_out.h>
43 #include <deal.II/numerics/derivative_approximation.h>
44
45 // And this again is C++:
46 #include <array>
47 #include <iostream>
48 #include <fstream>
49
50 // The last step is as in all previous programs:
51 namespace Step30
52 {
53 using namespace dealii;
54
55 // @sect3{Equation data}
56 //
57 // The classes describing equation data and the actual assembly of
58 // individual terms are almost entirely copied from step-12. We will comment
59 // on differences.
60 template <int dim>
61 class RHS : public Function<dim>
62 {
63 public:
value_list(const std::vector<Point<dim>> & points,std::vector<double> & values,const unsigned int=0) const64 virtual void value_list(const std::vector<Point<dim>> &points,
65 std::vector<double> & values,
66 const unsigned int /*component*/ = 0) const override
67 {
68 (void)points;
69 Assert(values.size() == points.size(),
70 ExcDimensionMismatch(values.size(), points.size()));
71
72 std::fill(values.begin(), values.end(), 0.);
73 }
74 };
75
76
77 template <int dim>
78 class BoundaryValues : public Function<dim>
79 {
80 public:
value_list(const std::vector<Point<dim>> & points,std::vector<double> & values,const unsigned int=0) const81 virtual void value_list(const std::vector<Point<dim>> &points,
82 std::vector<double> & values,
83 const unsigned int /*component*/ = 0) const override
84 {
85 Assert(values.size() == points.size(),
86 ExcDimensionMismatch(values.size(), points.size()));
87
88 for (unsigned int i = 0; i < values.size(); ++i)
89 {
90 if (points[i](0) < 0.5)
91 values[i] = 1.;
92 else
93 values[i] = 0.;
94 }
95 }
96 };
97
98
99 template <int dim>
100 class Beta
101 {
102 public:
103 // The flow field is chosen to be a quarter circle with counterclockwise
104 // flow direction and with the origin as midpoint for the right half of the
105 // domain with positive $x$ values, whereas the flow simply goes to the left
106 // in the left part of the domain at a velocity that matches the one coming
107 // in from the right. In the circular part the magnitude of the flow
108 // velocity is proportional to the distance from the origin. This is a
109 // difference to step-12, where the magnitude was 1 everywhere. the new
110 // definition leads to a linear variation of $\beta$ along each given face
111 // of a cell. On the other hand, the solution $u(x,y)$ is exactly the same
112 // as before.
value_list(const std::vector<Point<dim>> & points,std::vector<Point<dim>> & values) const113 void value_list(const std::vector<Point<dim>> &points,
114 std::vector<Point<dim>> & values) const
115 {
116 Assert(values.size() == points.size(),
117 ExcDimensionMismatch(values.size(), points.size()));
118
119 for (unsigned int i = 0; i < points.size(); ++i)
120 {
121 if (points[i](0) > 0)
122 {
123 values[i](0) = -points[i](1);
124 values[i](1) = points[i](0);
125 }
126 else
127 {
128 values[i] = Point<dim>();
129 values[i](0) = -points[i](1);
130 }
131 }
132 }
133 };
134
135
136
137 // @sect3{Class: DGTransportEquation}
138 //
139 // This declaration of this class is utterly unaffected by our current
140 // changes.
141 template <int dim>
142 class DGTransportEquation
143 {
144 public:
145 DGTransportEquation();
146
147 void assemble_cell_term(const FEValues<dim> &fe_v,
148 FullMatrix<double> & ui_vi_matrix,
149 Vector<double> & cell_vector) const;
150
151 void assemble_boundary_term(const FEFaceValues<dim> &fe_v,
152 FullMatrix<double> & ui_vi_matrix,
153 Vector<double> & cell_vector) const;
154
155 void assemble_face_term(const FEFaceValuesBase<dim> &fe_v,
156 const FEFaceValuesBase<dim> &fe_v_neighbor,
157 FullMatrix<double> & ui_vi_matrix,
158 FullMatrix<double> & ue_vi_matrix,
159 FullMatrix<double> & ui_ve_matrix,
160 FullMatrix<double> & ue_ve_matrix) const;
161
162 private:
163 const Beta<dim> beta_function;
164 const RHS<dim> rhs_function;
165 const BoundaryValues<dim> boundary_function;
166 };
167
168
169
170 // Likewise, the constructor of the class as well as the functions
171 // assembling the terms corresponding to cell interiors and boundary faces
172 // are unchanged from before. The function that assembles face terms between
173 // cells also did not change because all it does is operate on two objects
174 // of type FEFaceValuesBase (which is the base class of both FEFaceValues
175 // and FESubfaceValues). Where these objects come from, i.e. how they are
176 // initialized, is of no concern to this function: it simply assumes that
177 // the quadrature points on faces or subfaces represented by the two objects
178 // correspond to the same points in physical space.
179 template <int dim>
DGTransportEquation()180 DGTransportEquation<dim>::DGTransportEquation()
181 : beta_function()
182 , rhs_function()
183 , boundary_function()
184 {}
185
186
187
188 template <int dim>
assemble_cell_term(const FEValues<dim> & fe_v,FullMatrix<double> & ui_vi_matrix,Vector<double> & cell_vector) const189 void DGTransportEquation<dim>::assemble_cell_term(
190 const FEValues<dim> &fe_v,
191 FullMatrix<double> & ui_vi_matrix,
192 Vector<double> & cell_vector) const
193 {
194 const std::vector<double> &JxW = fe_v.get_JxW_values();
195
196 std::vector<Point<dim>> beta(fe_v.n_quadrature_points);
197 std::vector<double> rhs(fe_v.n_quadrature_points);
198
199 beta_function.value_list(fe_v.get_quadrature_points(), beta);
200 rhs_function.value_list(fe_v.get_quadrature_points(), rhs);
201
202 for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
203 for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
204 {
205 for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
206 ui_vi_matrix(i, j) -= beta[point] * fe_v.shape_grad(i, point) *
207 fe_v.shape_value(j, point) * JxW[point];
208
209 cell_vector(i) +=
210 rhs[point] * fe_v.shape_value(i, point) * JxW[point];
211 }
212 }
213
214
215 template <int dim>
assemble_boundary_term(const FEFaceValues<dim> & fe_v,FullMatrix<double> & ui_vi_matrix,Vector<double> & cell_vector) const216 void DGTransportEquation<dim>::assemble_boundary_term(
217 const FEFaceValues<dim> &fe_v,
218 FullMatrix<double> & ui_vi_matrix,
219 Vector<double> & cell_vector) const
220 {
221 const std::vector<double> & JxW = fe_v.get_JxW_values();
222 const std::vector<Tensor<1, dim>> &normals = fe_v.get_normal_vectors();
223
224 std::vector<Point<dim>> beta(fe_v.n_quadrature_points);
225 std::vector<double> g(fe_v.n_quadrature_points);
226
227 beta_function.value_list(fe_v.get_quadrature_points(), beta);
228 boundary_function.value_list(fe_v.get_quadrature_points(), g);
229
230 for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
231 {
232 const double beta_n = beta[point] * normals[point];
233 if (beta_n > 0)
234 for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
235 for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
236 ui_vi_matrix(i, j) += beta_n * fe_v.shape_value(j, point) *
237 fe_v.shape_value(i, point) * JxW[point];
238 else
239 for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
240 cell_vector(i) -=
241 beta_n * g[point] * fe_v.shape_value(i, point) * JxW[point];
242 }
243 }
244
245
246 template <int dim>
assemble_face_term(const FEFaceValuesBase<dim> & fe_v,const FEFaceValuesBase<dim> & fe_v_neighbor,FullMatrix<double> & ui_vi_matrix,FullMatrix<double> & ue_vi_matrix,FullMatrix<double> & ui_ve_matrix,FullMatrix<double> & ue_ve_matrix) const247 void DGTransportEquation<dim>::assemble_face_term(
248 const FEFaceValuesBase<dim> &fe_v,
249 const FEFaceValuesBase<dim> &fe_v_neighbor,
250 FullMatrix<double> & ui_vi_matrix,
251 FullMatrix<double> & ue_vi_matrix,
252 FullMatrix<double> & ui_ve_matrix,
253 FullMatrix<double> & ue_ve_matrix) const
254 {
255 const std::vector<double> & JxW = fe_v.get_JxW_values();
256 const std::vector<Tensor<1, dim>> &normals = fe_v.get_normal_vectors();
257
258 std::vector<Point<dim>> beta(fe_v.n_quadrature_points);
259
260 beta_function.value_list(fe_v.get_quadrature_points(), beta);
261
262 for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
263 {
264 const double beta_n = beta[point] * normals[point];
265 if (beta_n > 0)
266 {
267 for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
268 for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
269 ui_vi_matrix(i, j) += beta_n * fe_v.shape_value(j, point) *
270 fe_v.shape_value(i, point) * JxW[point];
271
272 for (unsigned int k = 0; k < fe_v_neighbor.dofs_per_cell; ++k)
273 for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
274 ui_ve_matrix(k, j) -= beta_n * fe_v.shape_value(j, point) *
275 fe_v_neighbor.shape_value(k, point) *
276 JxW[point];
277 }
278 else
279 {
280 for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
281 for (unsigned int l = 0; l < fe_v_neighbor.dofs_per_cell; ++l)
282 ue_vi_matrix(i, l) += beta_n *
283 fe_v_neighbor.shape_value(l, point) *
284 fe_v.shape_value(i, point) * JxW[point];
285
286 for (unsigned int k = 0; k < fe_v_neighbor.dofs_per_cell; ++k)
287 for (unsigned int l = 0; l < fe_v_neighbor.dofs_per_cell; ++l)
288 ue_ve_matrix(k, l) -=
289 beta_n * fe_v_neighbor.shape_value(l, point) *
290 fe_v_neighbor.shape_value(k, point) * JxW[point];
291 }
292 }
293 }
294
295
296 // @sect3{Class: DGMethod}
297 //
298 // This declaration is much like that of step-12. However, we introduce a
299 // new routine (set_anisotropic_flags) and modify another one (refine_grid).
300 template <int dim>
301 class DGMethod
302 {
303 public:
304 DGMethod(const bool anisotropic);
305
306 void run();
307
308 private:
309 void setup_system();
310 void assemble_system();
311 void solve(Vector<double> &solution);
312 void refine_grid();
313 void set_anisotropic_flags();
314 void output_results(const unsigned int cycle) const;
315
316 Triangulation<dim> triangulation;
317 const MappingQ1<dim> mapping;
318 // Again we want to use DG elements of degree 1 (but this is only
319 // specified in the constructor). If you want to use a DG method of a
320 // different degree replace 1 in the constructor by the new degree.
321 const unsigned int degree;
322 FE_DGQ<dim> fe;
323 DoFHandler<dim> dof_handler;
324
325 SparsityPattern sparsity_pattern;
326 SparseMatrix<double> system_matrix;
327 // This is new, the threshold value used in the evaluation of the
328 // anisotropic jump indicator explained in the introduction. Its value is
329 // set to 3.0 in the constructor, but it can easily be changed to a
330 // different value greater than 1.
331 const double anisotropic_threshold_ratio;
332 // This is a bool flag indicating whether anisotropic refinement shall be
333 // used or not. It is set by the constructor, which takes an argument of
334 // the same name.
335 const bool anisotropic;
336
337 const QGauss<dim> quadrature;
338 const QGauss<dim - 1> face_quadrature;
339
340 Vector<double> solution2;
341 Vector<double> right_hand_side;
342
343 const DGTransportEquation<dim> dg;
344 };
345
346
347 template <int dim>
DGMethod(const bool anisotropic)348 DGMethod<dim>::DGMethod(const bool anisotropic)
349 : mapping()
350 ,
351 // Change here for DG methods of different degrees.
352 degree(1)
353 , fe(degree)
354 , dof_handler(triangulation)
355 , anisotropic_threshold_ratio(3.)
356 , anisotropic(anisotropic)
357 ,
358 // As beta is a linear function, we can choose the degree of the
359 // quadrature for which the resulting integration is correct. Thus, we
360 // choose to use <code>degree+1</code> Gauss points, which enables us to
361 // integrate exactly polynomials of degree <code>2*degree+1</code>, enough
362 // for all the integrals we will perform in this program.
363 quadrature(degree + 1)
364 , face_quadrature(degree + 1)
365 , dg()
366 {}
367
368
369
370 template <int dim>
setup_system()371 void DGMethod<dim>::setup_system()
372 {
373 dof_handler.distribute_dofs(fe);
374 sparsity_pattern.reinit(dof_handler.n_dofs(),
375 dof_handler.n_dofs(),
376 (GeometryInfo<dim>::faces_per_cell *
377 GeometryInfo<dim>::max_children_per_face +
378 1) *
379 fe.n_dofs_per_cell());
380
381 DoFTools::make_flux_sparsity_pattern(dof_handler, sparsity_pattern);
382
383 sparsity_pattern.compress();
384
385 system_matrix.reinit(sparsity_pattern);
386
387 solution2.reinit(dof_handler.n_dofs());
388 right_hand_side.reinit(dof_handler.n_dofs());
389 }
390
391
392 // @sect4{Function: assemble_system}
393 //
394 // We proceed with the <code>assemble_system</code> function that implements
395 // the DG discretization. This function does the same thing as the
396 // <code>assemble_system</code> function from step-12 (but without
397 // MeshWorker). The four cases considered for the neighbor-relations of a
398 // cell are the same as the isotropic case, namely a) cell is at the
399 // boundary, b) there are finer neighboring cells, c) the neighbor is
400 // neither coarser nor finer and d) the neighbor is coarser. However, the
401 // way in which we decide upon which case we have are modified in the way
402 // described in the introduction.
403 template <int dim>
assemble_system()404 void DGMethod<dim>::assemble_system()
405 {
406 const unsigned int dofs_per_cell = dof_handler.get_fe().n_dofs_per_cell();
407 std::vector<types::global_dof_index> dofs(dofs_per_cell);
408 std::vector<types::global_dof_index> dofs_neighbor(dofs_per_cell);
409
410 const UpdateFlags update_flags = update_values | update_gradients |
411 update_quadrature_points |
412 update_JxW_values;
413
414 const UpdateFlags face_update_flags =
415 update_values | update_quadrature_points | update_JxW_values |
416 update_normal_vectors;
417
418 const UpdateFlags neighbor_face_update_flags = update_values;
419
420 FEValues<dim> fe_v(mapping, fe, quadrature, update_flags);
421 FEFaceValues<dim> fe_v_face(mapping,
422 fe,
423 face_quadrature,
424 face_update_flags);
425 FESubfaceValues<dim> fe_v_subface(mapping,
426 fe,
427 face_quadrature,
428 face_update_flags);
429 FEFaceValues<dim> fe_v_face_neighbor(mapping,
430 fe,
431 face_quadrature,
432 neighbor_face_update_flags);
433
434
435 FullMatrix<double> ui_vi_matrix(dofs_per_cell, dofs_per_cell);
436 FullMatrix<double> ue_vi_matrix(dofs_per_cell, dofs_per_cell);
437
438 FullMatrix<double> ui_ve_matrix(dofs_per_cell, dofs_per_cell);
439 FullMatrix<double> ue_ve_matrix(dofs_per_cell, dofs_per_cell);
440
441 Vector<double> cell_vector(dofs_per_cell);
442
443 for (const auto &cell : dof_handler.active_cell_iterators())
444 {
445 ui_vi_matrix = 0;
446 cell_vector = 0;
447
448 fe_v.reinit(cell);
449
450 dg.assemble_cell_term(fe_v, ui_vi_matrix, cell_vector);
451
452 cell->get_dof_indices(dofs);
453
454 for (const auto face_no : cell->face_indices())
455 {
456 const auto face = cell->face(face_no);
457
458 // Case (a): The face is at the boundary.
459 if (face->at_boundary())
460 {
461 fe_v_face.reinit(cell, face_no);
462
463 dg.assemble_boundary_term(fe_v_face, ui_vi_matrix, cell_vector);
464 }
465 else
466 {
467 Assert(cell->neighbor(face_no).state() == IteratorState::valid,
468 ExcInternalError());
469 const auto neighbor = cell->neighbor(face_no);
470
471 // Case (b): This is an internal face and the neighbor
472 // is refined (which we can test by asking whether the
473 // face of the current cell has children). In this
474 // case, we will need to integrate over the
475 // "sub-faces", i.e., the children of the face of the
476 // current cell.
477 //
478 // (There is a slightly confusing corner case: If we
479 // are in 1d -- where admittedly the current program
480 // and its demonstration of anisotropic refinement is
481 // not particularly relevant -- then the faces between
482 // cells are always the same: they are just
483 // vertices. In other words, in 1d, we do not want to
484 // treat faces between cells of different level
485 // differently. The condition `face->has_children()`
486 // we check here ensures this: in 1d, this function
487 // always returns `false`, and consequently in 1d we
488 // will not ever go into this `if` branch. But we will
489 // have to come back to this corner case below in case
490 // (c).)
491 if (face->has_children())
492 {
493 // We need to know, which of the neighbors faces points in
494 // the direction of our cell. Using the @p
495 // neighbor_face_no function we get this information for
496 // both coarser and non-coarser neighbors.
497 const unsigned int neighbor2 =
498 cell->neighbor_face_no(face_no);
499
500 // Now we loop over all subfaces, i.e. the children and
501 // possibly grandchildren of the current face.
502 for (unsigned int subface_no = 0;
503 subface_no < face->number_of_children();
504 ++subface_no)
505 {
506 // To get the cell behind the current subface we can
507 // use the @p neighbor_child_on_subface function. it
508 // takes care of all the complicated situations of
509 // anisotropic refinement and non-standard faces.
510 const auto neighbor_child =
511 cell->neighbor_child_on_subface(face_no, subface_no);
512 Assert(!neighbor_child->has_children(),
513 ExcInternalError());
514
515 // The remaining part of this case is unchanged.
516 ue_vi_matrix = 0;
517 ui_ve_matrix = 0;
518 ue_ve_matrix = 0;
519
520 fe_v_subface.reinit(cell, face_no, subface_no);
521 fe_v_face_neighbor.reinit(neighbor_child, neighbor2);
522
523 dg.assemble_face_term(fe_v_subface,
524 fe_v_face_neighbor,
525 ui_vi_matrix,
526 ue_vi_matrix,
527 ui_ve_matrix,
528 ue_ve_matrix);
529
530 neighbor_child->get_dof_indices(dofs_neighbor);
531
532 for (unsigned int i = 0; i < dofs_per_cell; ++i)
533 for (unsigned int j = 0; j < dofs_per_cell; ++j)
534 {
535 system_matrix.add(dofs[i],
536 dofs_neighbor[j],
537 ue_vi_matrix(i, j));
538 system_matrix.add(dofs_neighbor[i],
539 dofs[j],
540 ui_ve_matrix(i, j));
541 system_matrix.add(dofs_neighbor[i],
542 dofs_neighbor[j],
543 ue_ve_matrix(i, j));
544 }
545 }
546 }
547 else
548 {
549 // Case (c). We get here if this is an internal
550 // face and if the neighbor is not further refined
551 // (or, as mentioned above, we are in 1d in which
552 // case we get here for every internal face). We
553 // then need to decide whether we want to
554 // integrate over the current face. If the
555 // neighbor is in fact coarser, then we ignore the
556 // face and instead handle it when we visit the
557 // neighboring cell and look at the current face
558 // (except in 1d, where as mentioned above this is
559 // not happening):
560 if (dim > 1 && cell->neighbor_is_coarser(face_no))
561 continue;
562
563 // On the other hand, if the neighbor is more
564 // refined, then we have already handled the face
565 // in case (b) above (except in 1d). So for 2d and
566 // 3d, we just have to decide whether we want to
567 // handle a face between cells at the same level
568 // from the current side or from the neighboring
569 // side. We do this by introducing a tie-breaker:
570 // We'll just take the cell with the smaller index
571 // (within the current refinement level). In 1d,
572 // we take either the coarser cell, or if they are
573 // on the same level, the one with the smaller
574 // index within that level. This leads to a
575 // complicated condition that, hopefully, makes
576 // sense given the description above:
577 if (((dim > 1) && (cell->index() < neighbor->index())) ||
578 ((dim == 1) && ((cell->level() < neighbor->level()) ||
579 ((cell->level() == neighbor->level()) &&
580 (cell->index() < neighbor->index())))))
581 {
582 // Here we know, that the neighbor is not coarser so we
583 // can use the usual @p neighbor_of_neighbor
584 // function. However, we could also use the more
585 // general @p neighbor_face_no function.
586 const unsigned int neighbor2 =
587 cell->neighbor_of_neighbor(face_no);
588
589 ue_vi_matrix = 0;
590 ui_ve_matrix = 0;
591 ue_ve_matrix = 0;
592
593 fe_v_face.reinit(cell, face_no);
594 fe_v_face_neighbor.reinit(neighbor, neighbor2);
595
596 dg.assemble_face_term(fe_v_face,
597 fe_v_face_neighbor,
598 ui_vi_matrix,
599 ue_vi_matrix,
600 ui_ve_matrix,
601 ue_ve_matrix);
602
603 neighbor->get_dof_indices(dofs_neighbor);
604
605 for (unsigned int i = 0; i < dofs_per_cell; ++i)
606 for (unsigned int j = 0; j < dofs_per_cell; ++j)
607 {
608 system_matrix.add(dofs[i],
609 dofs_neighbor[j],
610 ue_vi_matrix(i, j));
611 system_matrix.add(dofs_neighbor[i],
612 dofs[j],
613 ui_ve_matrix(i, j));
614 system_matrix.add(dofs_neighbor[i],
615 dofs_neighbor[j],
616 ue_ve_matrix(i, j));
617 }
618 }
619
620 // We do not need to consider a case (d), as those
621 // faces are treated 'from the other side within
622 // case (b).
623 }
624 }
625 }
626
627 for (unsigned int i = 0; i < dofs_per_cell; ++i)
628 for (unsigned int j = 0; j < dofs_per_cell; ++j)
629 system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i, j));
630
631 for (unsigned int i = 0; i < dofs_per_cell; ++i)
632 right_hand_side(dofs[i]) += cell_vector(i);
633 }
634 }
635
636
637 // @sect3{Solver}
638 //
639 // For this simple problem we use the simple Richardson iteration again. The
640 // solver is completely unaffected by our anisotropic changes.
641 template <int dim>
solve(Vector<double> & solution)642 void DGMethod<dim>::solve(Vector<double> &solution)
643 {
644 SolverControl solver_control(1000, 1e-12, false, false);
645 SolverRichardson<Vector<double>> solver(solver_control);
646
647 PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
648
649 preconditioner.initialize(system_matrix, fe.n_dofs_per_cell());
650
651 solver.solve(system_matrix, solution, right_hand_side, preconditioner);
652 }
653
654
655 // @sect3{Refinement}
656 //
657 // We refine the grid according to the same simple refinement criterion used
658 // in step-12, namely an approximation to the gradient of the solution.
659 template <int dim>
refine_grid()660 void DGMethod<dim>::refine_grid()
661 {
662 Vector<float> gradient_indicator(triangulation.n_active_cells());
663
664 // We approximate the gradient,
665 DerivativeApproximation::approximate_gradient(mapping,
666 dof_handler,
667 solution2,
668 gradient_indicator);
669
670 // and scale it to obtain an error indicator.
671 for (const auto &cell : triangulation.active_cell_iterators())
672 gradient_indicator[cell->active_cell_index()] *=
673 std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
674 // Then we use this indicator to flag the 30 percent of the cells with
675 // highest error indicator to be refined.
676 GridRefinement::refine_and_coarsen_fixed_number(triangulation,
677 gradient_indicator,
678 0.3,
679 0.1);
680 // Now the refinement flags are set for those cells with a large error
681 // indicator. If nothing is done to change this, those cells will be
682 // refined isotropically. If the @p anisotropic flag given to this
683 // function is set, we now call the set_anisotropic_flags() function,
684 // which uses the jump indicator to reset some of the refinement flags to
685 // anisotropic refinement.
686 if (anisotropic)
687 set_anisotropic_flags();
688 // Now execute the refinement considering anisotropic as well as isotropic
689 // refinement flags.
690 triangulation.execute_coarsening_and_refinement();
691 }
692
693 // Once an error indicator has been evaluated and the cells with largest
694 // error are flagged for refinement we want to loop over the flagged cells
695 // again to decide whether they need isotropic refinement or whether
696 // anisotropic refinement is more appropriate. This is the anisotropic jump
697 // indicator explained in the introduction.
698 template <int dim>
set_anisotropic_flags()699 void DGMethod<dim>::set_anisotropic_flags()
700 {
701 // We want to evaluate the jump over faces of the flagged cells, so we
702 // need some objects to evaluate values of the solution on faces.
703 UpdateFlags face_update_flags =
704 UpdateFlags(update_values | update_JxW_values);
705
706 FEFaceValues<dim> fe_v_face(mapping,
707 fe,
708 face_quadrature,
709 face_update_flags);
710 FESubfaceValues<dim> fe_v_subface(mapping,
711 fe,
712 face_quadrature,
713 face_update_flags);
714 FEFaceValues<dim> fe_v_face_neighbor(mapping,
715 fe,
716 face_quadrature,
717 update_values);
718
719 // Now we need to loop over all active cells.
720 for (const auto &cell : dof_handler.active_cell_iterators())
721 // We only need to consider cells which are flagged for refinement.
722 if (cell->refine_flag_set())
723 {
724 Point<dim> jump;
725 Point<dim> area;
726
727 for (const auto face_no : cell->face_indices())
728 {
729 const auto face = cell->face(face_no);
730
731 if (!face->at_boundary())
732 {
733 Assert(cell->neighbor(face_no).state() ==
734 IteratorState::valid,
735 ExcInternalError());
736 const auto neighbor = cell->neighbor(face_no);
737
738 std::vector<double> u(fe_v_face.n_quadrature_points);
739 std::vector<double> u_neighbor(fe_v_face.n_quadrature_points);
740
741 // The four cases of different neighbor relations seen in
742 // the assembly routines are repeated much in the same way
743 // here.
744 if (face->has_children())
745 {
746 // The neighbor is refined. First we store the
747 // information, which of the neighbor's faces points in
748 // the direction of our current cell. This property is
749 // inherited to the children.
750 unsigned int neighbor2 = cell->neighbor_face_no(face_no);
751 // Now we loop over all subfaces,
752 for (unsigned int subface_no = 0;
753 subface_no < face->number_of_children();
754 ++subface_no)
755 {
756 // get an iterator pointing to the cell behind the
757 // present subface...
758 const auto neighbor_child =
759 cell->neighbor_child_on_subface(face_no,
760 subface_no);
761 Assert(!neighbor_child->has_children(),
762 ExcInternalError());
763 // ... and reinit the respective FEFaceValues and
764 // FESubFaceValues objects.
765 fe_v_subface.reinit(cell, face_no, subface_no);
766 fe_v_face_neighbor.reinit(neighbor_child, neighbor2);
767 // We obtain the function values
768 fe_v_subface.get_function_values(solution2, u);
769 fe_v_face_neighbor.get_function_values(solution2,
770 u_neighbor);
771 // as well as the quadrature weights, multiplied by
772 // the Jacobian determinant.
773 const std::vector<double> &JxW =
774 fe_v_subface.get_JxW_values();
775 // Now we loop over all quadrature points
776 for (unsigned int x = 0;
777 x < fe_v_subface.n_quadrature_points;
778 ++x)
779 {
780 // and integrate the absolute value of the jump
781 // of the solution, i.e. the absolute value of
782 // the difference between the function value
783 // seen from the current cell and the
784 // neighboring cell, respectively. We know, that
785 // the first two faces are orthogonal to the
786 // first coordinate direction on the unit cell,
787 // the second two faces are orthogonal to the
788 // second coordinate direction and so on, so we
789 // accumulate these values into vectors with
790 // <code>dim</code> components.
791 jump[face_no / 2] +=
792 std::abs(u[x] - u_neighbor[x]) * JxW[x];
793 // We also sum up the scaled weights to obtain
794 // the measure of the face.
795 area[face_no / 2] += JxW[x];
796 }
797 }
798 }
799 else
800 {
801 if (!cell->neighbor_is_coarser(face_no))
802 {
803 // Our current cell and the neighbor have the same
804 // refinement along the face under
805 // consideration. Apart from that, we do much the
806 // same as with one of the subcells in the above
807 // case.
808 unsigned int neighbor2 =
809 cell->neighbor_of_neighbor(face_no);
810
811 fe_v_face.reinit(cell, face_no);
812 fe_v_face_neighbor.reinit(neighbor, neighbor2);
813
814 fe_v_face.get_function_values(solution2, u);
815 fe_v_face_neighbor.get_function_values(solution2,
816 u_neighbor);
817
818 const std::vector<double> &JxW =
819 fe_v_face.get_JxW_values();
820
821 for (unsigned int x = 0;
822 x < fe_v_face.n_quadrature_points;
823 ++x)
824 {
825 jump[face_no / 2] +=
826 std::abs(u[x] - u_neighbor[x]) * JxW[x];
827 area[face_no / 2] += JxW[x];
828 }
829 }
830 else // i.e. neighbor is coarser than cell
831 {
832 // Now the neighbor is actually coarser. This case
833 // is new, in that it did not occur in the assembly
834 // routine. Here, we have to consider it, but this
835 // is not overly complicated. We simply use the @p
836 // neighbor_of_coarser_neighbor function, which
837 // again takes care of anisotropic refinement and
838 // non-standard face orientation by itself.
839 std::pair<unsigned int, unsigned int>
840 neighbor_face_subface =
841 cell->neighbor_of_coarser_neighbor(face_no);
842 Assert(neighbor_face_subface.first < cell->n_faces(),
843 ExcInternalError());
844 Assert(neighbor_face_subface.second <
845 neighbor->face(neighbor_face_subface.first)
846 ->number_of_children(),
847 ExcInternalError());
848 Assert(neighbor->neighbor_child_on_subface(
849 neighbor_face_subface.first,
850 neighbor_face_subface.second) == cell,
851 ExcInternalError());
852
853 fe_v_face.reinit(cell, face_no);
854 fe_v_subface.reinit(neighbor,
855 neighbor_face_subface.first,
856 neighbor_face_subface.second);
857
858 fe_v_face.get_function_values(solution2, u);
859 fe_v_subface.get_function_values(solution2,
860 u_neighbor);
861
862 const std::vector<double> &JxW =
863 fe_v_face.get_JxW_values();
864
865 for (unsigned int x = 0;
866 x < fe_v_face.n_quadrature_points;
867 ++x)
868 {
869 jump[face_no / 2] +=
870 std::abs(u[x] - u_neighbor[x]) * JxW[x];
871 area[face_no / 2] += JxW[x];
872 }
873 }
874 }
875 }
876 }
877 // Now we analyze the size of the mean jumps, which we get dividing
878 // the jumps by the measure of the respective faces.
879 std::array<double, dim> average_jumps;
880 double sum_of_average_jumps = 0.;
881 for (unsigned int i = 0; i < dim; ++i)
882 {
883 average_jumps[i] = jump(i) / area(i);
884 sum_of_average_jumps += average_jumps[i];
885 }
886
887 // Now we loop over the <code>dim</code> coordinate directions of
888 // the unit cell and compare the average jump over the faces
889 // orthogonal to that direction with the average jumps over faces
890 // orthogonal to the remaining direction(s). If the first is larger
891 // than the latter by a given factor, we refine only along hat
892 // axis. Otherwise we leave the refinement flag unchanged, resulting
893 // in isotropic refinement.
894 for (unsigned int i = 0; i < dim; ++i)
895 if (average_jumps[i] > anisotropic_threshold_ratio *
896 (sum_of_average_jumps - average_jumps[i]))
897 cell->set_refine_flag(RefinementCase<dim>::cut_axis(i));
898 }
899 }
900
901 // @sect3{The Rest}
902 //
903 // The remaining part of the program very much follows the scheme of
904 // previous tutorial programs. We output the mesh in VTU format (just
905 // as we did in step-1, for example), and the visualization output
906 // in VTU format as we almost always do.
907 template <int dim>
output_results(const unsigned int cycle) const908 void DGMethod<dim>::output_results(const unsigned int cycle) const
909 {
910 std::string refine_type;
911 if (anisotropic)
912 refine_type = ".aniso";
913 else
914 refine_type = ".iso";
915
916 {
917 const std::string filename =
918 "grid-" + std::to_string(cycle) + refine_type + ".svg";
919 std::cout << " Writing grid to <" << filename << ">..." << std::endl;
920 std::ofstream svg_output(filename);
921
922 GridOut grid_out;
923 grid_out.write_svg(triangulation, svg_output);
924 }
925
926 {
927 const std::string filename =
928 "sol-" + std::to_string(cycle) + refine_type + ".vtu";
929 std::cout << " Writing solution to <" << filename << ">..."
930 << std::endl;
931 std::ofstream gnuplot_output(filename);
932
933 DataOut<dim> data_out;
934 data_out.attach_dof_handler(dof_handler);
935 data_out.add_data_vector(solution2, "u");
936
937 data_out.build_patches(degree);
938
939 data_out.write_vtu(gnuplot_output);
940 }
941 }
942
943
944
945 template <int dim>
run()946 void DGMethod<dim>::run()
947 {
948 for (unsigned int cycle = 0; cycle < 6; ++cycle)
949 {
950 std::cout << "Cycle " << cycle << ':' << std::endl;
951
952 if (cycle == 0)
953 {
954 // Create the rectangular domain.
955 Point<dim> p1, p2;
956 p1(0) = 0;
957 p1(0) = -1;
958 for (unsigned int i = 0; i < dim; ++i)
959 p2(i) = 1.;
960 // Adjust the number of cells in different directions to obtain
961 // completely isotropic cells for the original mesh.
962 std::vector<unsigned int> repetitions(dim, 1);
963 repetitions[0] = 2;
964 GridGenerator::subdivided_hyper_rectangle(triangulation,
965 repetitions,
966 p1,
967 p2);
968
969 triangulation.refine_global(5 - dim);
970 }
971 else
972 refine_grid();
973
974
975 std::cout << " Number of active cells: "
976 << triangulation.n_active_cells() << std::endl;
977
978 setup_system();
979
980 std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
981 << std::endl;
982
983 Timer assemble_timer;
984 assemble_system();
985 std::cout << " Time of assemble_system: " << assemble_timer.cpu_time()
986 << std::endl;
987 solve(solution2);
988
989 output_results(cycle);
990
991 std::cout << std::endl;
992 }
993 }
994 } // namespace Step30
995
996
997
main()998 int main()
999 {
1000 try
1001 {
1002 using namespace Step30;
1003
1004 // If you want to run the program in 3D, simply change the following
1005 // line to <code>const unsigned int dim = 3;</code>.
1006 const unsigned int dim = 2;
1007
1008 {
1009 // First, we perform a run with isotropic refinement.
1010 std::cout << "Performing a " << dim
1011 << "D run with isotropic refinement..." << std::endl
1012 << "------------------------------------------------"
1013 << std::endl;
1014 DGMethod<dim> dgmethod_iso(false);
1015 dgmethod_iso.run();
1016 }
1017
1018 {
1019 // Now we do a second run, this time with anisotropic refinement.
1020 std::cout << std::endl
1021 << "Performing a " << dim
1022 << "D run with anisotropic refinement..." << std::endl
1023 << "--------------------------------------------------"
1024 << std::endl;
1025 DGMethod<dim> dgmethod_aniso(true);
1026 dgmethod_aniso.run();
1027 }
1028 }
1029 catch (std::exception &exc)
1030 {
1031 std::cerr << std::endl
1032 << std::endl
1033 << "----------------------------------------------------"
1034 << std::endl;
1035 std::cerr << "Exception on processing: " << std::endl
1036 << exc.what() << std::endl
1037 << "Aborting!" << std::endl
1038 << "----------------------------------------------------"
1039 << std::endl;
1040 return 1;
1041 }
1042 catch (...)
1043 {
1044 std::cerr << std::endl
1045 << std::endl
1046 << "----------------------------------------------------"
1047 << std::endl;
1048 std::cerr << "Unknown exception!" << std::endl
1049 << "Aborting!" << std::endl
1050 << "----------------------------------------------------"
1051 << std::endl;
1052 return 1;
1053 };
1054
1055 return 0;
1056 }
1057