--- /dev/null
+New: The step-12 tutorial program has been changed to use
+FEInterfaceValues. The old version of step-12 is still available as step-12b.
+<br>
+(Timo Heister, 2019/08/27)
+<br>
+
+<i> Note: A variant called step-12b of this tutorial exists, that uses
+MeshWorker and LocalIntegrators instead of assembling matrices using
+FEInterfaceValues as it is done in this tutorial.
+</i>
+
<a name="Intro"></a>
-<h1>An example of the MeshWorker framework with an advection problem</h1>
+<h1>An example of an advection problem using the Discountinuous Galerkin method</h1>
<h3>Overview</h3>
-This example is devoted to the MeshWorker framework and the <em>discontinuous
-Galerkin method</em>, or in short: DG method. It includes the following topics.
+This example is devoted to the <em>discontinuous
+Galerkin method</em>, or in short, the DG method. It includes the following topics.
<ol>
<li> Discretization of the linear advection equation with the DG method.
- <li> Assembling of the system matrix using the MeshWorker::loop().
+ <li> Assembling of jump terms and other expressions on the interface between cells using FEInterfaceValues.
+ <li> Assembling of the system matrix using the MeshWorker::mesh_loop().
</ol>
The particular concern of this program are the loops of DG methods. These turn
out to be especially complex, primarily because for the face terms, we have to
distinguish the cases of boundary, regular interior faces and interior faces
-with hanging nodes, respectively. The MeshWorker framework implements the
-standard loop over all cells and faces in MeshWorker::loop() and takes care of
-distinguishing between all the different faces.
-
-There are two things left to do if you use MeshWorker: first, you need
-to write the local integrators for your problem. Second, you select
-classes from the MeshWorker namespace and combine them to achieve your
-goal.
-
+with hanging nodes, respectively. The MeshWorker::mesh_loop() handles the
+complexity on iterating over cells and faces and allows specifying "workers"
+for the different cell and face terms. The integration of face terms itself,
+including on faces with hanging nodes, is done using the FEInterfaceValues
+class.
<h3>Problem</h3>
terms $\bigl(u_h,{\mathbf \beta}\cdot\nabla v_h\bigr)_T$, the internal fluxes
$\bigl<u_h^-, \beta\cdot[v_h\mathbf n]\bigr>_{F}$, and the boundary terms $\bigl<u_h,
v_h \beta\cdot \mathbf n]\bigr>_{\Gamma_+}$ and $\bigl<g, \beta\cdot\mathbf n
-v_h\bigr>_{\Gamma_-}$. The summation of all those is done by MeshWorker::integration_loop().
+v_h\bigr>_{\Gamma_-}$. The summation of all those is done by MeshWorker::mesh_loop().
*
* Author: Guido Kanschat, Texas A&M University, 2009
+ * Timo Heister, Clemson University, 2019
*/
// classes at all: they are passed to <code>DoFHandler</code> and
// <code>FEValues</code> objects, and that is about it.
#include <deal.II/fe/fe_dgq.h>
+// This header is needed for FEInterfaceValues to compute integrals on
+// interfaces
+#include <deal.II/fe/fe_interface_values.h>
// We are going to use the simplest possible solver, called Richardson
// iteration, that represents a simple defect correction. This, in combination
// with a block SSOR preconditioner (defined in precondition_block.h), that
// We are going to use gradients as refinement indicator.
#include <deal.II/numerics/derivative_approximation.h>
-// Here come the new include files for using the MeshWorker framework. The first
-// contains the class MeshWorker::DoFInfo, which provides local integrators with
-// a mapping between local and global degrees of freedom. It stores the results
-// of local integrals as well in its base class MeshWorker::LocalResults.
-// In the second of these files, we find an object of type
-// MeshWorker::IntegrationInfo, which is mostly a wrapper around a group of
-// FEValues objects. The file <tt>meshworker/simple.h</tt> contains classes
-// assembling locally integrated data into a global system containing only a
-// single matrix. Finally, we will need the file that runs the loop over all
-// mesh cells and faces.
-#include <deal.II/meshworker/dof_info.h>
-#include <deal.II/meshworker/integration_info.h>
-#include <deal.II/meshworker/simple.h>
-#include <deal.II/meshworker/loop.h>
+// Here come the new include files for using the mesh_loop from the MeshWorker
+// framework
+#include <deal.II/meshworker/mesh_loop.h>
// Like in all programs, we finish this section by including the needed C++
// headers and declaring we want to use objects in the dealii namespace without
private:
void setup_system();
void assemble_system();
- void solve(Vector<double> &solution);
+ void solve();
void refine_grid();
void output_results(const unsigned int cycle) const;
FE_DGQ<dim> fe;
DoFHandler<dim> dof_handler;
+ AffineConstraints<> constraints;
+
// The next four members represent the linear system to be solved.
// <code>system_matrix</code> and <code>right_hand_side</code> are generated
// by <code>assemble_system()</code>, the <code>solution</code> is computed
Vector<double> solution;
Vector<double> right_hand_side;
-
- // Finally, we have to provide functions that assemble the cell, boundary,
- // and inner face terms. Within the MeshWorker framework, the loop over all
- // cells and much of the setup of operations will be done outside this
- // class, so all we have to provide are these three operations. They will
- // then work on intermediate objects for which first, we here define
- // alias to the info objects handed to the local integration functions
- // in order to make our life easier below.
- using DoFInfo = MeshWorker::DoFInfo<dim>;
- using CellInfo = MeshWorker::IntegrationInfo<dim>;
-
- // The following three functions are then the ones that get called inside
- // the generic loop over all cells and faces. They are the ones doing the
- // actual integration.
- //
- // In our code below, these functions do not access member variables of the
- // current class, so we can mark them as <code>static</code> and simply pass
- // pointers to these functions to the MeshWorker framework. If, however,
- // these functions would want to access member variables (or needed
- // additional arguments beyond the ones specified below), we could use the
- // facilities of lambda functions to provide the
- // MeshWorker framework with objects that act as if they had the required
- // number and types of arguments, but have in fact other arguments already
- // bound.
- static void integrate_cell_term(DoFInfo &dinfo, CellInfo &info);
- static void integrate_boundary_term(DoFInfo &dinfo, CellInfo &info);
- static void integrate_face_term(DoFInfo & dinfo1,
- DoFInfo & dinfo2,
- CellInfo &info1,
- CellInfo &info2);
};
template <int dim>
void AdvectionProblem<dim>::assemble_system()
{
- // This is the magic object, which knows everything about the data
- // structures and local integration. This is the object doing the work in
- // the function MeshWorker::loop(), which is implicitly called by
- // MeshWorker::integration_loop() below. After the functions to which we
- // provide pointers did the local integration, the
- // MeshWorker::Assembler::SystemSimple object distributes these into the
- // global sparse matrix and the right hand side vector.
- MeshWorker::IntegrationInfoBox<dim> info_box;
-
- // First, we initialize the quadrature formulae and the update flags in the
- // worker base class. For quadrature, we play safe and use a QGauss formula
- // with number of points one higher than the polynomial degree used. Since
- // the quadratures for cells, boundary and interior faces can be selected
- // independently, we have to hand over this value three times.
- const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
- info_box.initialize_gauss_quadrature(n_gauss_points,
- n_gauss_points,
- n_gauss_points);
-
- // These are the types of values we need for integrating our system. They
- // are added to the flags used on cells, boundary and interior faces, as
- // well as interior neighbor faces, which is forced by the four @p true
- // values.
- info_box.initialize_update_flags();
- UpdateFlags update_flags =
- update_quadrature_points | update_values | update_gradients;
- info_box.add_update_flags(update_flags, true, true, true, true);
-
- // After preparing all data in <tt>info_box</tt>, we initialize the FEValues
- // objects in there.
- info_box.initialize(fe, mapping);
-
- // The object created so far helps us do the local integration on each cell
- // and face. Now, we need an object which receives the integrated (local)
- // data and forwards them to the assembler.
- MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
- // Now, we have to create the assembler object and tell it, where to put the
- // local data. These will be our system matrix and the right hand side.
- MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double>>
- assembler;
- assembler.initialize(system_matrix, right_hand_side);
-
- // Finally, the integration loop over all active cells (determined by the
- // first argument, which is an active iterator).
- //
- // As noted in the discussion when declaring the local integration functions
- // in the class declaration, the arguments expected by the assembling
- // integrator class are not actually function pointers. Rather, they are
- // objects that can be called like functions with a certain number of
- // arguments. Consequently, we could also pass objects with appropriate
- // operator() implementations here, or lambda functions if the local
- // integrators were, for example, non-static member functions.
- MeshWorker::loop<dim,
- dim,
- MeshWorker::DoFInfo<dim>,
- MeshWorker::IntegrationInfoBox<dim>>(
- dof_handler.begin_active(),
- dof_handler.end(),
- dof_info,
- info_box,
- &AdvectionProblem<dim>::integrate_cell_term,
- &AdvectionProblem<dim>::integrate_boundary_term,
- &AdvectionProblem<dim>::integrate_face_term,
- assembler);
- }
-
-
- // @sect4{The local integrators}
+ typedef decltype(dof_handler.begin_active()) Iterator;
+ BoundaryValues<dim> boundary_function;
+
+ auto cell_worker = [&](const Iterator & cell,
+ ScratchData<dim> &scratch_data,
+ CopyData & copy_data) {
+ const unsigned int n_dofs = scratch_data.fe_values.get_fe().dofs_per_cell;
+ copy_data.reinit(cell, n_dofs);
+ scratch_data.fe_values.reinit(cell);
+
+ const auto &q_points = scratch_data.fe_values.get_quadrature_points();
+
+ const FEValues<dim> & fe_v = scratch_data.fe_values;
+ const std::vector<double> &JxW = fe_v.get_JxW_values();
+
+ // We solve a homogeneous equation, thus no right hand side shows up in
+ // the cell term. What's left is integrating the matrix entries.
+ for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
+ {
+ auto beta_q = beta(q_points[point]);
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ for (unsigned int j = 0; j < n_dofs; ++j)
+ {
+ copy_data.cell_matrix(i, j) +=
+ -beta_q * fe_v.shape_grad(i, point) *
+ fe_v.shape_value(j, point) * JxW[point];
+ }
+ }
+ };
+
+ /*
+ auto boundary_worker1 = [&](const Iterator & cell,
+ const unsigned int &face_no,
+ ScratchData<dim> & scratch_data,
+ CopyData & copy_data) {
+ FEFacetValues<dim> &fe_facet = scratch_data.fe_facet_values;
+ fe_facet.reinit(cell, face_no);
+
+ const auto &q_points =
+ fe_facet.get_fe_values().get_quadrature_points();
+
+ const unsigned int n_facet_dofs = fe_facet.n_facet_dofs();
+ const std::vector<double> &JxW =
+ fe_facet.get_fe_values().get_JxW_values();
+ const std::vector<Tensor<1, dim>> &normals =
+ fe_facet.get_fe_values().get_normal_vectors();
+
+ std::vector<double> g(q_points.size());
+ boundary_function.value_list(q_points, g);
+
+ for (unsigned int point = 0; point < q_points.size(); ++point)
+ {
+ const double beta_n = beta(q_points[point]) * normals[point];
+
+ if (beta_n > 0)
+ {
+ for (unsigned int i = 0; i < n_facet_dofs / 2; ++i) // TODO:
+ ugly! for (unsigned int j = 0; j < n_facet_dofs / 2; ++j)
+ copy_data.cell_matrix(i, j) +=
+ beta_n * fe_facet.scalar().jump(j, point) *
+ fe_facet.scalar().choose(true, i, point) * JxW[point];
+ }
+ else if (0)
+ for (unsigned int i = 0; i < n_facet_dofs / 2; ++i) // TODO:
+ ugly copy_data.cell_rhs(i) -= beta_n * g[point] *
+ fe_facet.scalar().jump(i, point) *
+ JxW[point];
+ }
+ };*/
+
+ auto boundary_worker = [&](const Iterator & cell,
+ const unsigned int &face_no,
+ ScratchData<dim> & scratch_data,
+ CopyData & copy_data) {
+ scratch_data.fe_interface_values.reinit(cell, face_no);
+ const FEFaceValuesBase<dim> &fe_face =
+ scratch_data.fe_interface_values.get_fe_values();
+
+ const auto &q_points = fe_face.get_quadrature_points();
+
+ const unsigned int n_facet_dofs = fe_face.get_fe().n_dofs_per_cell();
+ const std::vector<double> & JxW = fe_face.get_JxW_values();
+ const std::vector<Tensor<1, dim>> &normals = fe_face.get_normal_vectors();
+
+ std::vector<double> g(q_points.size());
+ boundary_function.value_list(q_points, g);
+
+ for (unsigned int point = 0; point < q_points.size(); ++point)
+ {
+ const double beta_n = beta(q_points[point]) * normals[point];
+
+ if (beta_n > 0)
+ {
+ for (unsigned int i = 0; i < n_facet_dofs; ++i)
+ for (unsigned int j = 0; j < n_facet_dofs; ++j)
+ copy_data.cell_matrix(i, j) +=
+ beta_n * fe_face.shape_value(j, point) *
+ fe_face.shape_value(i, point) * JxW[point];
+ }
+ else
+ for (unsigned int i = 0; i < n_facet_dofs; ++i)
+ copy_data.cell_rhs(i) -=
+ beta_n * g[point] * fe_face.shape_value(i, point) * JxW[point];
+ }
+ };
+
+ auto face_worker = [&](const Iterator & cell,
+ const unsigned int &f,
+ const unsigned int &sf,
+ const Iterator & ncell,
+ const unsigned int &nf,
+ const unsigned int &nsf,
+ ScratchData<dim> & scratch_data,
+ CopyData & copy_data) {
+ FEInterfaceValues<dim> &fe_facet = scratch_data.fe_interface_values;
+ fe_facet.reinit(cell, f, sf, ncell, nf, nsf);
+ const auto &q_points = fe_facet.get_fe_values().get_quadrature_points();
+
+ copy_data.face_data.emplace_back();
+ CopyDataFace ©_data_face = copy_data.face_data.back();
+
+ const unsigned int n_dofs = fe_facet.n_interface_dofs();
+ copy_data_face.joint_dof_indices = fe_facet.get_interface_dof_indices();
+
+ copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
+
+ const std::vector<double> & JxW = fe_facet.get_JxW_values();
+ const std::vector<Tensor<1, dim>> &normals =
+ fe_facet.get_normal_vectors();
+
+ // u- * (beta*n)[v]
+ for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+ {
+ const double beta_n = beta(q_points[qpoint]) * normals[qpoint];
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ for (unsigned int j = 0; j < n_dofs; ++j)
+ copy_data_face.cell_matrix(i, j) +=
+ fe_facet.choose(beta_n > 0, j, qpoint) // u^-
+ * beta_n // (beta*n)
+ * fe_facet.jump(i, qpoint) // [v]
+ * JxW[qpoint]; // dx
+ }
+ };
+
+ auto copier = [&](const CopyData &c) {
+ constraints.distribute_local_to_global(c.cell_matrix,
+ c.cell_rhs,
+ c.local_dof_indices,
+ system_matrix,
+ right_hand_side);
+
+ for (auto &cdf : c.face_data)
+ {
+ constraints.distribute_local_to_global(cdf.cell_matrix,
+ cdf.joint_dof_indices,
+ system_matrix);
+ }
+ };
- // These are the functions given to the MeshWorker::integration_loop() called
- // just above. They compute the local contributions to the system matrix and
- // right hand side on cells and faces.
- template <int dim>
- void AdvectionProblem<dim>::integrate_cell_term(DoFInfo & dinfo,
- CellInfo &info)
- {
- // First, let us retrieve some of the objects used here from @p info. Note
- // that these objects can handle much more complex structures, thus the
- // access here looks more complicated than might seem necessary.
- const FEValuesBase<dim> & fe_values = info.fe_values();
- FullMatrix<double> & local_matrix = dinfo.matrix(0).matrix;
- const std::vector<double> &JxW = fe_values.get_JxW_values();
-
- // With these objects, we continue local integration like always. First, we
- // loop over the quadrature points and compute the advection vector in the
- // current point.
- for (unsigned int point = 0; point < fe_values.n_quadrature_points; ++point)
- {
- const Tensor<1, dim> beta_at_q_point =
- beta(fe_values.quadrature_point(point));
-
- // We solve a homogeneous equation, thus no right hand side shows up in
- // the cell term. What's left is integrating the matrix entries.
- for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
- for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
- local_matrix(i, j) += -beta_at_q_point * //
- fe_values.shape_grad(i, point) * //
- fe_values.shape_value(j, point) * //
- JxW[point];
- }
- }
-
- // Now the same for the boundary terms. Note that now we use FEValuesBase, the
- // base class for both FEFaceValues and FESubfaceValues, in order to get
- // access to normal vectors.
- template <int dim>
- void AdvectionProblem<dim>::integrate_boundary_term(DoFInfo & dinfo,
- CellInfo &info)
- {
- const FEValuesBase<dim> &fe_face_values = info.fe_values();
- FullMatrix<double> & local_matrix = dinfo.matrix(0).matrix;
- Vector<double> & local_vector = dinfo.vector(0).block(0);
-
- const std::vector<double> & JxW = fe_face_values.get_JxW_values();
- const std::vector<Tensor<1, dim>> &normals =
- fe_face_values.get_normal_vectors();
-
- std::vector<double> g(fe_face_values.n_quadrature_points);
-
- static BoundaryValues<dim> boundary_function;
- boundary_function.value_list(fe_face_values.get_quadrature_points(), g);
-
- for (unsigned int point = 0; point < fe_face_values.n_quadrature_points;
- ++point)
- {
- const double beta_dot_n =
- beta(fe_face_values.quadrature_point(point)) * normals[point];
- if (beta_dot_n > 0)
- for (unsigned int i = 0; i < fe_face_values.dofs_per_cell; ++i)
- for (unsigned int j = 0; j < fe_face_values.dofs_per_cell; ++j)
- local_matrix(i, j) += beta_dot_n * //
- fe_face_values.shape_value(j, point) * //
- fe_face_values.shape_value(i, point) * //
- JxW[point];
- else
- for (unsigned int i = 0; i < fe_face_values.dofs_per_cell; ++i)
- local_vector(i) += -beta_dot_n * //
- g[point] * //
- fe_face_values.shape_value(i, point) * //
- JxW[point];
- }
- }
+ const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
- // Finally, the interior face terms. The difference here is that we receive
- // two info objects, one for each cell adjacent to the face and we assemble
- // four matrices, one for each cell and two for coupling back and forth.
- template <int dim>
- void AdvectionProblem<dim>::integrate_face_term(DoFInfo & dinfo1,
- DoFInfo & dinfo2,
- CellInfo &info1,
- CellInfo &info2)
- {
- // For quadrature points, weights, etc., we use the FEValuesBase object of
- // the first argument.
- const FEValuesBase<dim> &fe_face_values = info1.fe_values();
- const unsigned int dofs_per_cell = fe_face_values.dofs_per_cell;
-
- // For additional shape functions, we have to ask the neighbors
- // FEValuesBase.
- const FEValuesBase<dim> &fe_face_values_neighbor = info2.fe_values();
- const unsigned int neighbor_dofs_per_cell =
- fe_face_values_neighbor.dofs_per_cell;
-
- // Then we get references to the four local matrices. The letters u and v
- // refer to trial and test functions, respectively. The %numbers indicate
- // the cells provided by info1 and info2. By convention, the two matrices
- // in each info object refer to the test functions on the respective cell.
- // The first matrix contains the interior couplings of that cell, while the
- // second contains the couplings between cells.
- FullMatrix<double> &u1_v1_matrix = dinfo1.matrix(0, false).matrix;
- FullMatrix<double> &u2_v1_matrix = dinfo1.matrix(0, true).matrix;
- FullMatrix<double> &u1_v2_matrix = dinfo2.matrix(0, true).matrix;
- FullMatrix<double> &u2_v2_matrix = dinfo2.matrix(0, false).matrix;
-
- // Here, following the previous functions, we would have the local right
- // hand side vectors. Fortunately, the interface terms only involve the
- // solution and the right hand side does not receive any contributions.
-
- const std::vector<double> & JxW = fe_face_values.get_JxW_values();
- const std::vector<Tensor<1, dim>> &normals =
- fe_face_values.get_normal_vectors();
-
- for (unsigned int point = 0; point < fe_face_values.n_quadrature_points;
- ++point)
- {
- const double beta_dot_n =
- beta(fe_face_values.quadrature_point(point)) * normals[point];
- if (beta_dot_n > 0)
- {
- // This term we've already seen:
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
- u1_v1_matrix(i, j) += beta_dot_n * //
- fe_face_values.shape_value(j, point) * //
- fe_face_values.shape_value(i, point) * //
- JxW[point];
-
- // We additionally assemble the term $(\beta\cdot n u,\hat
- // v)_{\partial \kappa_+}$,
- for (unsigned int k = 0; k < neighbor_dofs_per_cell; ++k)
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
- u1_v2_matrix(k, j) +=
- -beta_dot_n * //
- fe_face_values.shape_value(j, point) * //
- fe_face_values_neighbor.shape_value(k, point) * //
- JxW[point];
- }
- else
- {
- // This one we've already seen, too:
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- for (unsigned int l = 0; l < neighbor_dofs_per_cell; ++l)
- u2_v1_matrix(i, l) +=
- beta_dot_n * //
- fe_face_values_neighbor.shape_value(l, point) * //
- fe_face_values.shape_value(i, point) * //
- JxW[point];
-
- // And this is another new one: $(\beta\cdot n \hat u,\hat
- // v)_{\partial \kappa_-}$:
- for (unsigned int k = 0; k < neighbor_dofs_per_cell; ++k)
- for (unsigned int l = 0; l < neighbor_dofs_per_cell; ++l)
- u2_v2_matrix(k, l) +=
- -beta_dot_n * //
- fe_face_values_neighbor.shape_value(l, point) * //
- fe_face_values_neighbor.shape_value(k, point) * //
- JxW[point];
- }
- }
+ ScratchData<dim> scratch_data(mapping, fe, n_gauss_points);
+ CopyData copy_data;
+ MeshWorker::mesh_loop(dof_handler.begin_active(),
+ dof_handler.end(),
+ cell_worker,
+ copier,
+ scratch_data,
+ copy_data,
+ MeshWorker::assemble_own_cells |
+ MeshWorker::assemble_boundary_faces |
+ MeshWorker::assemble_own_interior_faces_once,
+ boundary_worker,
+ face_worker);
}
-
// @sect3{All the rest}
//
// For this simple problem we use the simplest possible solver, called
// then a block Gauss-Seidel preconditioner (see the PreconditionBlockSOR
// class with relaxation=1) does a much better job.
template <int dim>
- void AdvectionProblem<dim>::solve(Vector<double> &solution)
+ void AdvectionProblem<dim>::solve()
{
SolverControl solver_control(1000, 1e-12);
SolverRichardson<> solver(solver_control);
template <int dim>
void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
{
- // First write the grid in eps format.
- {
- const std::string filename = "grid-" + std::to_string(cycle) + ".eps";
- deallog << "Writing grid to <" << filename << ">" << std::endl;
- std::ofstream eps_output(filename);
-
- GridOut grid_out;
- grid_out.write_eps(triangulation, eps_output);
- }
-
- // Then output the solution in gnuplot format.
- {
- const std::string filename = "sol-" + std::to_string(cycle) + ".gnuplot";
- deallog << "Writing solution to <" << filename << ">" << std::endl;
- std::ofstream gnuplot_output(filename);
+ const std::string filename = "solution-" + std::to_string(cycle) + ".vtk";
+ deallog << "Writing solution to <" << filename << ">" << std::endl;
+ std::ofstream output(filename);
- DataOut<dim> data_out;
- data_out.attach_dof_handler(dof_handler);
- data_out.add_data_vector(solution, "u");
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "u");
- data_out.build_patches();
+ data_out.build_patches();
- data_out.write_gnuplot(gnuplot_output);
- }
+ data_out.write_vtk(output);
}
if (cycle == 0)
{
GridGenerator::hyper_cube(triangulation);
-
triangulation.refine_global(3);
}
else
refine_grid();
-
deallog << "Number of active cells: "
<< triangulation.n_active_cells() << std::endl;
<< std::endl;
assemble_system();
- solve(solution);
+ solve();
output_results(cycle);
}
// well, and need not be commented on.
int main()
{
+ dealii::deallog.depth_console(10);
try
{
- dealii::deallog.depth_console(5);
-
Step12::AdvectionProblem<2> dgmethod;
dgmethod.run();
}
--- /dev/null
+##
+# CMake script for the step-12b tutorial program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-12b")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
+# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+SET(TARGET_SRC
+ ${TARGET}.cc
+ )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+
+FIND_PACKAGE(deal.II 9.2.0 QUIET
+ HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+ )
+IF(NOT ${deal.II_FOUND})
+ MESSAGE(FATAL_ERROR "\n"
+ "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+ "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+ "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+ )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
--- /dev/null
+step-16 step-7 step-39
--- /dev/null
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+This is a variant of step-16 with the only change that we are using the
+MeshWorker framework with the pre-made LocalIntegrator helper classes instead
+of assembling the face terms using FEInterfaceValues.
+
+The details of this framework on how it is used in practice will be explained
+as part of this tutorial program.
+
+<h3>The testcase</h3>
+
+The problem we solve here is the same as the one in step-12.
+
--- /dev/null
+techniques
--- /dev/null
+<h1>Results</h1>
+
+
+The output of this program is very similar to step-16 and we are not repeating the output here.
+
+We show the solutions on the initial mesh, the mesh after two
+and after five adaptive refinement steps.
+
+<img src="https://www.dealii.org/images/steps/developer/step-12.sol-0.png" alt="">
+<img src="https://www.dealii.org/images/steps/developer/step-12.sol-2.png" alt="">
+<img src="https://www.dealii.org/images/steps/developer/step-12.sol-5.png" alt="">
+
+
+Then we show the final grid (after 5 refinement steps) and the solution again,
+this time with a nicer 3d rendering (obtained using the DataOutBase::write_vtk
+function and the VTK-based VisIt visualization program) that better shows the
+sharpness of the jump on the refined mesh and the over- and undershoots of the
+solution along the interface:
+
+<img src="https://www.dealii.org/images/steps/developer/step-12.grid-5.png" alt="">
+<img src="https://www.dealii.org/images/steps/developer/step-12.3d-solution.png" alt="">
+
+
+And finally we show a plot of a 3d computation.
+
+<img src="https://www.dealii.org/images/steps/developer/step-12.sol-5-3d.png" alt="">
+
+
+<a name="dg-vs-cg"></a>
+<h3>Why use discontinuous elements</h3>
+
+In this program we have used discontinuous elements. It is a legitimate
+question to ask why not simply use the normal, continuous ones. Of course, to
+everyone with a background in numerical methods, the answer is obvious: the
+continuous Galerkin (cG) method is not stable for the transport equation,
+unless one specifically adds stabilization terms. The DG method, however,
+<i>is</i> stable. Illustrating this with the current program is not very
+difficult; in fact, only the following minor modifications are necessary:
+- Change the element to FE_Q instead of FE_DGQ.
+- Add handling of hanging node constraints in exactly the same way as step-6.
+- We need a different solver; the direct solver in step-29 is a convenient
+ choice.
+An experienced deal.II user will be able to do this in less than 10 minutes.
+
+While the 2d solution has been shown above, containing a number of small
+spikes at the interface that are, however, stable in height under mesh
+refinement, results look much different when using a continuous element:
+
+<table align="center">
+ <tr>
+ <td valign="top">
+ 0
+ </td>
+ <td valign="top">
+ <img src="https://www.dealii.org/images/steps/developer/step-12.cg.sol-0.png" alt="">
+ </td>
+ <td valign="top">
+ 1
+ </td>
+ <td valign="top">
+ <img src="https://www.dealii.org/images/steps/developer/step-12.cg.sol-1.png" alt="">
+ </td>
+ </tr>
+ <tr>
+ <td valign="top">
+ 2
+ </td>
+ <td valign="top">
+ <img src="https://www.dealii.org/images/steps/developer/step-12.cg.sol-2.png" alt="">
+ </td>
+ <td valign="top">
+ 3
+ </td>
+ <td valign="top">
+ <img src="https://www.dealii.org/images/steps/developer/step-12.cg.sol-3.png" alt="">
+ </td>
+ </tr>
+ <tr>
+ <td valign="top">
+ 4
+ </td>
+ <td valign="top">
+ <img src="https://www.dealii.org/images/steps/developer/step-12.cg.sol-4.png" alt="">
+ </td>
+ <td valign="top">
+ 5
+ </td>
+ <td valign="top">
+ <img src="https://www.dealii.org/images/steps/developer/step-12.cg.sol-5.png" alt="">
+ </td>
+ </tr>
+</table>
+
+In refinement iteration 5, the image can't be plotted in a reasonable way any
+more as a 3d plot. We thus show a color plot with a range of $[-1,2]$ (the
+solution values of the exact solution lie in $[0,1]$, of course). In any case,
+it is clear that the continuous Galerkin solution exhibits oscillatory
+behavior that gets worse and worse as the mesh is refined more and more.
+
+There are a number of strategies to stabilize the cG method, if one wants to
+use continuous elements for some reason. Discussing these methods is beyond
+the scope of this tutorial program; an interested reader could, for example,
+take a look at step-31.
+
+
+
+<a name="extensions"></a>
+<h3>Possibilities for extensions</h3>
+
+Given that the exact solution is known in this case, one interesting
+avenue for further extensions would be to confirm the order of
+convergence for this program. In the current case, the solution is
+non-smooth, and so we can not expect to get a particularly high order
+of convergence, even if we used higher order elements. But even if the
+solution <i>is</i> smooth, the equation is not elliptic and so it is not
+immediately clear that we should obtain a convergence order that
+equals that of the optimal interpolation estimates (i.e. for example
+that we would get $h^3$ convergence in the $L^2$ norm by using
+quadratic elements).
+
+In fact, for hyperbolic equations, theoretical predictions often
+indicate that the best one can hope for is an order one half below the
+interpolation estimate. For example, for the streamline diffusion
+method (an alternative method to the DG method used here to stabilize
+the solution of the transport equation), one can prove that for
+elements of degree $p$, the order of convergence is $p+\frac 12$ on
+arbitrary meshes. While the observed order is frequently $p+1$ on
+uniformly refined meshes, one can construct so-called Peterson meshes
+on which the worse theoretical bound is actually attained. This should
+be relatively simple to verify, for example using the
+VectorTools::integrate_difference function.
+
+A different direction is to observe that the solution of transport problems
+often has discontinuities and that therefore a mesh in which we <i>bisect</i>
+every cell in every coordinate direction may not be optimal. Rather, a better
+strategy would be to only cut cells in the direction parallel to the
+discontinuity. This is called <i>anisotropic mesh refinement</i> and is the
+subject of step-30.
--- /dev/null
+Discontinuous Galerkin for linear advection, MeshWorker version.
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2009 - 2018 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Author: Guido Kanschat, Texas A&M University, 2009
+ */
+
+
+// The first few files have already been covered in previous examples and will
+// thus not be further commented on:
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/fe/mapping_q1.h>
+// Here the discontinuous finite elements are defined. They are used in the same
+// way as all other finite elements, though -- as you have seen in previous
+// tutorial programs -- there isn't much user interaction with finite element
+// classes at all: they are passed to <code>DoFHandler</code> and
+// <code>FEValues</code> objects, and that is about it.
+#include <deal.II/fe/fe_dgq.h>
+// We are going to use the simplest possible solver, called Richardson
+// iteration, that represents a simple defect correction. This, in combination
+// with a block SSOR preconditioner (defined in precondition_block.h), that
+// uses the special block matrix structure of system matrices arising from DG
+// discretizations.
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/precondition_block.h>
+// We are going to use gradients as refinement indicator.
+#include <deal.II/numerics/derivative_approximation.h>
+
+// Here come the new include files for using the MeshWorker framework. The first
+// contains the class MeshWorker::DoFInfo, which provides local integrators with
+// a mapping between local and global degrees of freedom. It stores the results
+// of local integrals as well in its base class MeshWorker::LocalResults.
+// In the second of these files, we find an object of type
+// MeshWorker::IntegrationInfo, which is mostly a wrapper around a group of
+// FEValues objects. The file <tt>meshworker/simple.h</tt> contains classes
+// assembling locally integrated data into a global system containing only a
+// single matrix. Finally, we will need the file that runs the loop over all
+// mesh cells and faces.
+#include <deal.II/meshworker/dof_info.h>
+#include <deal.II/meshworker/integration_info.h>
+#include <deal.II/meshworker/simple.h>
+#include <deal.II/meshworker/loop.h>
+
+// Like in all programs, we finish this section by including the needed C++
+// headers and declaring we want to use objects in the dealii namespace without
+// prefix.
+#include <iostream>
+#include <fstream>
+
+
+namespace Step12
+{
+ using namespace dealii;
+
+ // @sect3{Equation data}
+ //
+ // First, we define a class describing the inhomogeneous boundary data. Since
+ // only its values are used, we implement value_list(), but leave all other
+ // functions of Function undefined.
+ template <int dim>
+ class BoundaryValues : public Function<dim>
+ {
+ public:
+ BoundaryValues() = default;
+ virtual void value_list(const std::vector<Point<dim>> &points,
+ std::vector<double> & values,
+ const unsigned int component = 0) const override;
+ };
+
+ // Given the flow direction, the inflow boundary of the unit square $[0,1]^2$
+ // are the right and the lower boundaries. We prescribe discontinuous boundary
+ // values 1 and 0 on the x-axis and value 0 on the right boundary. The values
+ // of this function on the outflow boundaries will not be used within the DG
+ // scheme.
+ template <int dim>
+ void BoundaryValues<dim>::value_list(const std::vector<Point<dim>> &points,
+ std::vector<double> & values,
+ const unsigned int component) const
+ {
+ (void)component;
+ AssertIndexRange(component, 1);
+ Assert(values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i = 0; i < values.size(); ++i)
+ {
+ if (points[i](0) < 0.5)
+ values[i] = 1.;
+ else
+ values[i] = 0.;
+ }
+ }
+
+
+ // Finally, a function that computes and returns the wind field
+ // $\beta=\beta(\mathbf x)$. As explained in the introduction, we will use a
+ // rotational field around the origin in 2d. In 3d, we simply leave the
+ // $z$-component unset (i.e., at zero), whereas the function can not be used
+ // in 1d in its current implementation:
+ template <int dim>
+ Tensor<1, dim> beta(const Point<dim> &p)
+ {
+ Assert(dim >= 2, ExcNotImplemented());
+
+ Point<dim> wind_field;
+ wind_field(0) = -p(1);
+ wind_field(1) = p(0);
+ wind_field /= wind_field.norm();
+
+ return wind_field;
+ }
+
+
+ // @sect3{The AdvectionProblem class}
+ //
+ // After this preparations, we proceed with the main class of this program,
+ // called AdvectionProblem. It is basically the main class of step-6. We do
+ // not have an AffineConstraints object, because there are no hanging node
+ // constraints in DG discretizations.
+
+ // Major differences will only come up in the implementation of the assemble
+ // functions, since here, we not only need to cover the flux integrals over
+ // faces, we also use the MeshWorker interface to simplify the loops
+ // involved.
+ template <int dim>
+ class AdvectionProblem
+ {
+ public:
+ AdvectionProblem();
+ void run();
+
+ private:
+ void setup_system();
+ void assemble_system();
+ void solve(Vector<double> &solution);
+ void refine_grid();
+ void output_results(const unsigned int cycle) const;
+
+ Triangulation<dim> triangulation;
+ const MappingQ1<dim> mapping;
+
+ // Furthermore we want to use DG elements of degree 1 (but this is only
+ // specified in the constructor). If you want to use a DG method of a
+ // different degree the whole program stays the same, only replace 1 in
+ // the constructor by the desired polynomial degree.
+ FE_DGQ<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ // The next four members represent the linear system to be solved.
+ // <code>system_matrix</code> and <code>right_hand_side</code> are generated
+ // by <code>assemble_system()</code>, the <code>solution</code> is computed
+ // in <code>solve()</code>. The <code>sparsity_pattern</code> is used to
+ // determine the location of nonzero elements in <code>system_matrix</code>.
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> right_hand_side;
+
+ // Finally, we have to provide functions that assemble the cell, boundary,
+ // and inner face terms. Within the MeshWorker framework, the loop over all
+ // cells and much of the setup of operations will be done outside this
+ // class, so all we have to provide are these three operations. They will
+ // then work on intermediate objects for which first, we here define
+ // alias to the info objects handed to the local integration functions
+ // in order to make our life easier below.
+ using DoFInfo = MeshWorker::DoFInfo<dim>;
+ using CellInfo = MeshWorker::IntegrationInfo<dim>;
+
+ // The following three functions are then the ones that get called inside
+ // the generic loop over all cells and faces. They are the ones doing the
+ // actual integration.
+ //
+ // In our code below, these functions do not access member variables of the
+ // current class, so we can mark them as <code>static</code> and simply pass
+ // pointers to these functions to the MeshWorker framework. If, however,
+ // these functions would want to access member variables (or needed
+ // additional arguments beyond the ones specified below), we could use the
+ // facilities of boost::bind (or std::bind, respectively) to provide the
+ // MeshWorker framework with objects that act as if they had the required
+ // number and types of arguments, but have in fact other arguments already
+ // bound.
+ static void integrate_cell_term(DoFInfo &dinfo, CellInfo &info);
+ static void integrate_boundary_term(DoFInfo &dinfo, CellInfo &info);
+ static void integrate_face_term(DoFInfo & dinfo1,
+ DoFInfo & dinfo2,
+ CellInfo &info1,
+ CellInfo &info2);
+ };
+
+
+ // We start with the constructor. The 1 in the constructor call of
+ // <code>fe</code> is the polynomial degree.
+ template <int dim>
+ AdvectionProblem<dim>::AdvectionProblem()
+ : mapping()
+ , fe(1)
+ , dof_handler(triangulation)
+ {}
+
+
+ template <int dim>
+ void AdvectionProblem<dim>::setup_system()
+ {
+ // In the function that sets up the usual finite element data structures, we
+ // first need to distribute the DoFs.
+ dof_handler.distribute_dofs(fe);
+
+ // We start by generating the sparsity pattern. To this end, we first fill
+ // an intermediate object of type DynamicSparsityPattern with the couplings
+ // appearing in the system. After building the pattern, this object is
+ // copied to <code>sparsity_pattern</code> and can be discarded.
+
+ // To build the sparsity pattern for DG discretizations, we can call the
+ // function analogue to DoFTools::make_sparsity_pattern, which is called
+ // DoFTools::make_flux_sparsity_pattern:
+ DynamicSparsityPattern dsp(dof_handler.n_dofs());
+ DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
+ sparsity_pattern.copy_from(dsp);
+
+ // Finally, we set up the structure of all components of the linear system.
+ system_matrix.reinit(sparsity_pattern);
+ solution.reinit(dof_handler.n_dofs());
+ right_hand_side.reinit(dof_handler.n_dofs());
+ }
+
+ // @sect4{The assemble_system function}
+
+ // Here we see the major difference to assembling by hand. Instead of writing
+ // loops over cells and faces, we leave all this to the MeshWorker framework.
+ // In order to do so, we just have to define local integration functions and
+ // use one of the classes in namespace MeshWorker::Assembler to build the
+ // global system.
+ template <int dim>
+ void AdvectionProblem<dim>::assemble_system()
+ {
+ // This is the magic object, which knows everything about the data
+ // structures and local integration. This is the object doing the work in
+ // the function MeshWorker::loop(), which is implicitly called by
+ // MeshWorker::integration_loop() below. After the functions to which we
+ // provide pointers did the local integration, the
+ // MeshWorker::Assembler::SystemSimple object distributes these into the
+ // global sparse matrix and the right hand side vector.
+ MeshWorker::IntegrationInfoBox<dim> info_box;
+
+ // First, we initialize the quadrature formulae and the update flags in the
+ // worker base class. For quadrature, we play safe and use a QGauss formula
+ // with number of points one higher than the polynomial degree used. Since
+ // the quadratures for cells, boundary and interior faces can be selected
+ // independently, we have to hand over this value three times.
+ const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
+ info_box.initialize_gauss_quadrature(n_gauss_points,
+ n_gauss_points,
+ n_gauss_points);
+
+ // These are the types of values we need for integrating our system. They
+ // are added to the flags used on cells, boundary and interior faces, as
+ // well as interior neighbor faces, which is forced by the four @p true
+ // values.
+ info_box.initialize_update_flags();
+ UpdateFlags update_flags =
+ update_quadrature_points | update_values | update_gradients;
+ info_box.add_update_flags(update_flags, true, true, true, true);
+
+ // After preparing all data in <tt>info_box</tt>, we initialize the FEValues
+ // objects in there.
+ info_box.initialize(fe, mapping);
+
+ // The object created so far helps us do the local integration on each cell
+ // and face. Now, we need an object which receives the integrated (local)
+ // data and forwards them to the assembler.
+ MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+
+ // Now, we have to create the assembler object and tell it, where to put the
+ // local data. These will be our system matrix and the right hand side.
+ MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double>>
+ assembler;
+ assembler.initialize(system_matrix, right_hand_side);
+
+ // Finally, the integration loop over all active cells (determined by the
+ // first argument, which is an active iterator).
+ //
+ // As noted in the discussion when declaring the local integration functions
+ // in the class declaration, the arguments expected by the assembling
+ // integrator class are not actually function pointers. Rather, they are
+ // objects that can be called like functions with a certain number of
+ // arguments. Consequently, we could also pass objects with appropriate
+ // operator() implementations here, or the result of std::bind if the local
+ // integrators were, for example, non-static member functions.
+ MeshWorker::loop<dim,
+ dim,
+ MeshWorker::DoFInfo<dim>,
+ MeshWorker::IntegrationInfoBox<dim>>(
+ dof_handler.begin_active(),
+ dof_handler.end(),
+ dof_info,
+ info_box,
+ &AdvectionProblem<dim>::integrate_cell_term,
+ &AdvectionProblem<dim>::integrate_boundary_term,
+ &AdvectionProblem<dim>::integrate_face_term,
+ assembler);
+ }
+
+
+ // @sect4{The local integrators}
+
+ // These are the functions given to the MeshWorker::integration_loop() called
+ // just above. They compute the local contributions to the system matrix and
+ // right hand side on cells and faces.
+ template <int dim>
+ void AdvectionProblem<dim>::integrate_cell_term(DoFInfo & dinfo,
+ CellInfo &info)
+ {
+ // First, let us retrieve some of the objects used here from @p info. Note
+ // that these objects can handle much more complex structures, thus the
+ // access here looks more complicated than might seem necessary.
+ const FEValuesBase<dim> & fe_values = info.fe_values();
+ FullMatrix<double> & local_matrix = dinfo.matrix(0).matrix;
+ const std::vector<double> &JxW = fe_values.get_JxW_values();
+
+ // With these objects, we continue local integration like always. First, we
+ // loop over the quadrature points and compute the advection vector in the
+ // current point.
+ for (unsigned int point = 0; point < fe_values.n_quadrature_points; ++point)
+ {
+ const Tensor<1, dim> beta_at_q_point =
+ beta(fe_values.quadrature_point(point));
+
+ // We solve a homogeneous equation, thus no right hand side shows up in
+ // the cell term. What's left is integrating the matrix entries.
+ for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
+ local_matrix(i, j) += -beta_at_q_point * //
+ fe_values.shape_grad(i, point) * //
+ fe_values.shape_value(j, point) * //
+ JxW[point];
+ }
+ }
+
+ // Now the same for the boundary terms. Note that now we use FEValuesBase, the
+ // base class for both FEFaceValues and FESubfaceValues, in order to get
+ // access to normal vectors.
+ template <int dim>
+ void AdvectionProblem<dim>::integrate_boundary_term(DoFInfo & dinfo,
+ CellInfo &info)
+ {
+ const FEValuesBase<dim> &fe_face_values = info.fe_values();
+ FullMatrix<double> & local_matrix = dinfo.matrix(0).matrix;
+ Vector<double> & local_vector = dinfo.vector(0).block(0);
+
+ const std::vector<double> & JxW = fe_face_values.get_JxW_values();
+ const std::vector<Tensor<1, dim>> &normals =
+ fe_face_values.get_normal_vectors();
+
+ std::vector<double> g(fe_face_values.n_quadrature_points);
+
+ static BoundaryValues<dim> boundary_function;
+ boundary_function.value_list(fe_face_values.get_quadrature_points(), g);
+
+ for (unsigned int point = 0; point < fe_face_values.n_quadrature_points;
+ ++point)
+ {
+ const double beta_dot_n =
+ beta(fe_face_values.quadrature_point(point)) * normals[point];
+ if (beta_dot_n > 0)
+ for (unsigned int i = 0; i < fe_face_values.dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < fe_face_values.dofs_per_cell; ++j)
+ local_matrix(i, j) += beta_dot_n * //
+ fe_face_values.shape_value(j, point) * //
+ fe_face_values.shape_value(i, point) * //
+ JxW[point];
+ else
+ for (unsigned int i = 0; i < fe_face_values.dofs_per_cell; ++i)
+ local_vector(i) += -beta_dot_n * //
+ g[point] * //
+ fe_face_values.shape_value(i, point) * //
+ JxW[point];
+ }
+ }
+
+ // Finally, the interior face terms. The difference here is that we receive
+ // two info objects, one for each cell adjacent to the face and we assemble
+ // four matrices, one for each cell and two for coupling back and forth.
+ template <int dim>
+ void AdvectionProblem<dim>::integrate_face_term(DoFInfo & dinfo1,
+ DoFInfo & dinfo2,
+ CellInfo &info1,
+ CellInfo &info2)
+ {
+ // For quadrature points, weights, etc., we use the FEValuesBase object of
+ // the first argument.
+ const FEValuesBase<dim> &fe_face_values = info1.fe_values();
+ const unsigned int dofs_per_cell = fe_face_values.dofs_per_cell;
+
+ // For additional shape functions, we have to ask the neighbors
+ // FEValuesBase.
+ const FEValuesBase<dim> &fe_face_values_neighbor = info2.fe_values();
+ const unsigned int neighbor_dofs_per_cell =
+ fe_face_values_neighbor.dofs_per_cell;
+
+ // Then we get references to the four local matrices. The letters u and v
+ // refer to trial and test functions, respectively. The %numbers indicate
+ // the cells provided by info1 and info2. By convention, the two matrices
+ // in each info object refer to the test functions on the respective cell.
+ // The first matrix contains the interior couplings of that cell, while the
+ // second contains the couplings between cells.
+ FullMatrix<double> &u1_v1_matrix = dinfo1.matrix(0, false).matrix;
+ FullMatrix<double> &u2_v1_matrix = dinfo1.matrix(0, true).matrix;
+ FullMatrix<double> &u1_v2_matrix = dinfo2.matrix(0, true).matrix;
+ FullMatrix<double> &u2_v2_matrix = dinfo2.matrix(0, false).matrix;
+
+ // Here, following the previous functions, we would have the local right
+ // hand side vectors. Fortunately, the interface terms only involve the
+ // solution and the right hand side does not receive any contributions.
+
+ const std::vector<double> & JxW = fe_face_values.get_JxW_values();
+ const std::vector<Tensor<1, dim>> &normals =
+ fe_face_values.get_normal_vectors();
+
+ for (unsigned int point = 0; point < fe_face_values.n_quadrature_points;
+ ++point)
+ {
+ const double beta_dot_n =
+ beta(fe_face_values.quadrature_point(point)) * normals[point];
+ if (beta_dot_n > 0)
+ {
+ // This term we've already seen:
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ u1_v1_matrix(i, j) += beta_dot_n * //
+ fe_face_values.shape_value(j, point) * //
+ fe_face_values.shape_value(i, point) * //
+ JxW[point];
+
+ // We additionally assemble the term $(\beta\cdot n u,\hat
+ // v)_{\partial \kappa_+}$,
+ for (unsigned int k = 0; k < neighbor_dofs_per_cell; ++k)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ u1_v2_matrix(k, j) +=
+ -beta_dot_n * //
+ fe_face_values.shape_value(j, point) * //
+ fe_face_values_neighbor.shape_value(k, point) * //
+ JxW[point];
+ }
+ else
+ {
+ // This one we've already seen, too:
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int l = 0; l < neighbor_dofs_per_cell; ++l)
+ u2_v1_matrix(i, l) +=
+ beta_dot_n * //
+ fe_face_values_neighbor.shape_value(l, point) * //
+ fe_face_values.shape_value(i, point) * //
+ JxW[point];
+
+ // And this is another new one: $(\beta\cdot n \hat u,\hat
+ // v)_{\partial \kappa_-}$:
+ for (unsigned int k = 0; k < neighbor_dofs_per_cell; ++k)
+ for (unsigned int l = 0; l < neighbor_dofs_per_cell; ++l)
+ u2_v2_matrix(k, l) +=
+ -beta_dot_n * //
+ fe_face_values_neighbor.shape_value(l, point) * //
+ fe_face_values_neighbor.shape_value(k, point) * //
+ JxW[point];
+ }
+ }
+ }
+
+
+ // @sect3{All the rest}
+ //
+ // For this simple problem we use the simplest possible solver, called
+ // Richardson iteration, that represents a simple defect correction. This, in
+ // combination with a block SSOR preconditioner, that uses the special block
+ // matrix structure of system matrices arising from DG discretizations. The
+ // size of these blocks are the number of DoFs per cell. Here, we use a SSOR
+ // preconditioning as we have not renumbered the DoFs according to the flow
+ // field. If the DoFs are renumbered in the downstream direction of the flow,
+ // then a block Gauss-Seidel preconditioner (see the PreconditionBlockSOR
+ // class with relaxation=1) does a much better job.
+ template <int dim>
+ void AdvectionProblem<dim>::solve(Vector<double> &solution)
+ {
+ SolverControl solver_control(1000, 1e-12);
+ SolverRichardson<> solver(solver_control);
+
+ // Here we create the preconditioner,
+ PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
+
+ // then assign the matrix to it and set the right block size:
+ preconditioner.initialize(system_matrix, fe.dofs_per_cell);
+
+ // After these preparations we are ready to start the linear solver.
+ solver.solve(system_matrix, solution, right_hand_side, preconditioner);
+ }
+
+
+ // We refine the grid according to a very simple refinement criterion, namely
+ // an approximation to the gradient of the solution. As here we consider the
+ // DG(1) method (i.e. we use piecewise bilinear shape functions) we could
+ // simply compute the gradients on each cell. But we do not want to base our
+ // refinement indicator on the gradients on each cell only, but want to base
+ // them also on jumps of the discontinuous solution function over faces
+ // between neighboring cells. The simplest way of doing that is to compute
+ // approximative gradients by difference quotients including the cell under
+ // consideration and its neighbors. This is done by the
+ // <code>DerivativeApproximation</code> class that computes the approximate
+ // gradients in a way similar to the <code>GradientEstimation</code> described
+ // in step-9 of this tutorial. In fact, the
+ // <code>DerivativeApproximation</code> class was developed following the
+ // <code>GradientEstimation</code> class of step-9. Relating to the discussion
+ // in step-9, here we consider $h^{1+d/2}|\nabla_h u_h|$. Furthermore we note
+ // that we do not consider approximate second derivatives because solutions to
+ // the linear advection equation are in general not in $H^2$ but only in $H^1$
+ // (or, to be more precise: in $H^1_\beta$, i.e., the space of functions whose
+ // derivatives in direction $\beta$ are square integrable).
+ template <int dim>
+ void AdvectionProblem<dim>::refine_grid()
+ {
+ // The <code>DerivativeApproximation</code> class computes the gradients to
+ // float precision. This is sufficient as they are approximate and serve as
+ // refinement indicators only.
+ Vector<float> gradient_indicator(triangulation.n_active_cells());
+
+ // Now the approximate gradients are computed
+ DerivativeApproximation::approximate_gradient(mapping,
+ dof_handler,
+ solution,
+ gradient_indicator);
+
+ // and they are cell-wise scaled by the factor $h^{1+d/2}$
+ unsigned int cell_no = 0;
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ gradient_indicator(cell_no++) *=
+ std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
+
+ // Finally they serve as refinement indicator.
+ GridRefinement::refine_and_coarsen_fixed_number(triangulation,
+ gradient_indicator,
+ 0.3,
+ 0.1);
+
+ triangulation.execute_coarsening_and_refinement();
+ }
+
+
+ // The output of this program consists of eps-files of the adaptively refined
+ // grids and the numerical solutions given in gnuplot format.
+ template <int dim>
+ void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
+ {
+ // First write the grid in eps format.
+ {
+ const std::string filename = "grid-" + std::to_string(cycle) + ".eps";
+ deallog << "Writing grid to <" << filename << ">" << std::endl;
+ std::ofstream eps_output(filename);
+
+ GridOut grid_out;
+ grid_out.write_eps(triangulation, eps_output);
+ }
+
+ // Then output the solution in gnuplot format.
+ {
+ const std::string filename = "sol-" + std::to_string(cycle) + ".gnuplot";
+ deallog << "Writing solution to <" << filename << ">" << std::endl;
+ std::ofstream gnuplot_output(filename);
+
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "u");
+
+ data_out.build_patches();
+
+ data_out.write_gnuplot(gnuplot_output);
+ }
+ }
+
+
+ // The following <code>run</code> function is similar to previous examples.
+ template <int dim>
+ void AdvectionProblem<dim>::run()
+ {
+ for (unsigned int cycle = 0; cycle < 6; ++cycle)
+ {
+ deallog << "Cycle " << cycle << std::endl;
+
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_cube(triangulation);
+
+ triangulation.refine_global(3);
+ }
+ else
+ refine_grid();
+
+
+ deallog << "Number of active cells: "
+ << triangulation.n_active_cells() << std::endl;
+
+ setup_system();
+
+ deallog << "Number of degrees of freedom: " << dof_handler.n_dofs()
+ << std::endl;
+
+ assemble_system();
+ solve(solution);
+
+ output_results(cycle);
+ }
+ }
+} // namespace Step12
+
+
+// The following <code>main</code> function is similar to previous examples as
+// well, and need not be commented on.
+int main()
+{
+ try
+ {
+ dealii::deallog.depth_console(5);
+
+ Step12::AdvectionProblem<2> dgmethod;
+ dgmethod.run();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl
+ << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl
+ << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+
+ return 0;
+}
projects / dealii.git / commitdiff