Table of contents | |
---|---|
This example is devoted to the MeshWorker framework and the discontinuous Galerkin method, or in short: DG method. It includes the following topics.
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.
The model problem solved in this example is the linear advection equation
subject to the boundary conditions
on the inflow part of the boundary
of the domain. Here,
denotes a vector field,
the (scalar) solution function,
a boundary value function,
the inflow part of the boundary of the domain and denotes the unit outward normal to the boundary
. This equation is the conservative version of the advection equation already considered in step-9 of this tutorial. In particular, we solve the advection equation on
with
representing a circular counterclockwise flow field, and
on
and
on
.
We apply the well-known upwind discontinuous Galerkin method. To this end, we introduce the mesh dependent bilinear form
Here, is the set of all active cells of the triangulation and
is the set of all active interior faces.
and
denote the L2-inner products on the cell
and a face
, respectively. The jump is defined as
, where the superscripts refer to the upwind ('+') and downwind ('-') values at the face.
In order to implement this bilinear form, we need to compute the cell terms , the internal fluxes
, and the boundary terms
and
. The summation of all those is done by MeshWorker::integration_loop().
The first few files have already been covered in previous examples and will thus not be further commented on:
#include <base/quadrature_lib.h> #include <base/function.h> #include <lac/vector.h> #include <lac/compressed_sparsity_pattern.h> #include <lac/sparse_matrix.h> #include <grid/tria.h> #include <grid/grid_generator.h> #include <grid/grid_out.h> #include <grid/grid_refinement.h> #include <grid/tria_accessor.h> #include <grid/tria_iterator.h> #include <fe/fe_values.h> #include <dofs/dof_handler.h> #include <dofs/dof_accessor.h> #include <dofs/dof_tools.h> #include <numerics/data_out.h> #include <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: the are passed to DoFHandler
and FEValues
objects, and that is about it.
#include <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 <lac/solver_richardson.h> #include <lac/precondition_block.h>
We are going to use gradients as refinement indicator.
#include <numerics/derivative_approximation.h>
Here come the new include files for using the MeshWorker framework:
#include <numerics/mesh_worker.h> #include <numerics/mesh_worker_info.h> #include <numerics/mesh_worker_assembler.h> #include <numerics/mesh_worker_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> using namespace dealii;
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 () {}; virtual void value_list (const std::vector<Point<dim> > &points, std::vector<double> &values, const unsigned int component=0) const; };
Given the flow direction, the inflow boundary of the unit square 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) const { 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.; } }
After this preparations, we proceed with the main class of this program, called Step12. It is basically the main class of step-6. We do not have a ConstraintMatrix, 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 Step12 { public: Step12 (); 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. system_matrix
and right_hand_side
are generated by assemble_system()
, the solution
is computed in solve()
. The sparsity_pattern
is used to determine the location of nonzero elements in system_matrix
.
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 typedefs to the info objects handed to the local integration functions in order to make our life easier below.
typedef MeshWorker::DoFInfo<dim> DoFInfo; typedef MeshWorker::IntegrationInfo<dim> CellInfo;
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 static
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 fe
is the polynomial degree.
template <int dim> Step12<dim>::Step12 () : fe (1), dof_handler (triangulation) {} template <int dim> void Step12<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 CompressedSparsityPattern with the couplings appearing in the system. After building the pattern, this object is copied to sparsity_pattern
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:
CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); DoFTools::make_flux_sparsity_pattern (dof_handler, c_sparsity); sparsity_pattern.copy_from(c_sparsity);
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()); }
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 Step12<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 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 info_box
, we initialize the FEValus 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::integration_loop<dim, dim>
(dof_handler.begin_active(), dof_handler.end(),
dof_info, info_box,
&Step12<dim>::integrate_cell_term,
&Step12<dim>::integrate_boundary_term,
&Step12<dim>::integrate_face_term,
assembler, true);
}
These functions are analogous to step-12 and differ only in the data structures. Instead of providing the local matrices explicitly in the argument list, they are part of the info object.
Note that here we still have the local integration loop inside the following functions. The program would be even shorter, if we used pre-made operators from the Operators namespace (which will be added soon).
template <int dim> void Step12<dim>::integrate_cell_term (DoFInfo& dinfo, CellInfo& info) {
First, let us retrieve some of the objects used here from 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_v = info.fe_values(); FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix; const std::vector<double> &JxW = fe_v.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_v.n_quadrature_points; ++point) { Point<dim> beta; beta(0) = -fe_v.quadrature_point(point)(1); beta(1) = fe_v.quadrature_point(point)(0); beta /= beta.norm();
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_v.dofs_per_cell; ++i) for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j) local_matrix(i,j) -= beta*fe_v.shape_grad(i,point)* fe_v.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 Step12<dim>::integrate_boundary_term (DoFInfo& dinfo, CellInfo& info) { const FEValuesBase<dim>& fe_v = 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_v.get_JxW_values (); const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors (); std::vector<double> g(fe_v.n_quadrature_points); static BoundaryValues<dim> boundary_function; boundary_function.value_list (fe_v.get_quadrature_points(), g); for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point) { Point<dim> beta; beta(0) = -fe_v.quadrature_point(point)(1); beta(1) = fe_v.quadrature_point(point)(0); beta /= beta.norm(); const double beta_n=beta * normals[point]; if (beta_n>0) for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j) local_matrix(i,j) += beta_n * fe_v.shape_value(j,point) * fe_v.shape_value(i,point) * JxW[point]; else for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) local_vector(i) -= beta_n * g[point] * fe_v.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 Step12<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_v = info1.fe_values();
For additional shape functions, we have to ask the neighbors FEValuesBase.
const FEValuesBase<dim>& fe_v_neighbor = info2.fe_values();
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_v.get_JxW_values (); const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors (); for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point) { Point<dim> beta; beta(0) = -fe_v.quadrature_point(point)(1); beta(1) = fe_v.quadrature_point(point)(0); beta /= beta.norm(); const double beta_n=beta * normals[point]; if (beta_n>0) {
This term we've already seen:
for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j) u1_v1_matrix(i,j) += beta_n * fe_v.shape_value(j,point) * fe_v.shape_value(i,point) * JxW[point];
We additionally assemble the term ,
for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k) for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j) u1_v2_matrix(k,j) -= beta_n * fe_v.shape_value(j,point) * fe_v_neighbor.shape_value(k,point) * JxW[point]; } else {
This one we've already seen, too:
for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l) u2_v1_matrix(i,l) += beta_n * fe_v_neighbor.shape_value(l,point) * fe_v.shape_value(i,point) * JxW[point];
And this is another new one: :
for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k) for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l) u2_v2_matrix(k,l) -= beta_n * fe_v_neighbor.shape_value(l,point) * fe_v_neighbor.shape_value(k,point) * JxW[point]; } } }
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 Step12<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 DerivativeApproximation
class that computes the approximate gradients in a way similar to the GradientEstimation
described in step-9 of this tutorial. In fact, the DerivativeApproximation
class was developed following the GradientEstimation
class of step-9. Relating to the discussion in step-9, here we consider . Furthermore we note that we do not consider approximate second derivatives because solutions to the linear advection equation are in general not in
but in
(to be more precise, in
) only.
template <int dim> void Step12<dim>::refine_grid () {
The DerivativeApproximation
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
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) 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. This was covered in previous examples and will not be further commented on.
template <int dim> void Step12<dim>::output_results (const unsigned int cycle) const {
Write the grid in eps format.
std::string filename = "grid-"; filename += ('0' + cycle); Assert (cycle < 10, ExcInternalError()); filename += ".eps"; deallog << "Writing grid to <" << filename << ">" << std::endl; std::ofstream eps_output (filename.c_str()); GridOut grid_out; grid_out.write_eps (triangulation, eps_output);
Output of the solution in gnuplot format.
filename = "sol-"; filename += ('0' + cycle); Assert (cycle < 10, ExcInternalError()); filename += ".gnuplot"; deallog << "Writing solution to <" << filename << ">" << std::endl; std::ofstream gnuplot_output (filename.c_str()); 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 run
function is similar to previous examples.
template <int dim> void Step12<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); } }
The following main
function is similar to previous examples as well, and need not be commented on.
int main () { try { Step12<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; }
The output of this program consist of the console output, the eps files including the grids, and the solutions given in gnuplot format.
DEAL::Cycle 0 DEAL::Number of active cells: 64 DEAL::Number of degrees of freedom: 256 DEAL:Richardson::Starting value 0.176777 DEAL:Richardson::Convergence step 4 value 3.33123e-17 DEAL::Writing grid to <grid-0.eps> DEAL::Writing solution to <sol-0.gnuplot> DEAL::Cycle 1 DEAL::Number of active cells: 112 DEAL::Number of degrees of freedom: 448 DEAL:Richardson::Starting value 0.153093 DEAL:Richardson::Convergence step 9 value 3.74479e-17 DEAL::Writing grid to <grid-1.eps> DEAL::Writing solution to <sol-1.gnuplot> DEAL::Cycle 2 DEAL::Number of active cells: 214 DEAL::Number of degrees of freedom: 856 DEAL:Richardson::Starting value 0.149870 DEAL:Richardson::Convergence step 16 value 1.41017e-14 DEAL::Writing grid to <grid-2.eps> DEAL::Writing solution to <sol-2.gnuplot> DEAL::Cycle 3 DEAL::Number of active cells: 415 DEAL::Number of degrees of freedom: 1660 DEAL:Richardson::Starting value 0.149053 DEAL:Richardson::Convergence step 26 value 4.92424e-15 DEAL::Writing grid to <grid-3.eps> DEAL::Writing solution to <sol-3.gnuplot> DEAL::Cycle 4 DEAL::Number of active cells: 796 DEAL::Number of degrees of freedom: 3184 DEAL:Richardson::Starting value 0.148848 DEAL:Richardson::Convergence step 44 value 5.80787e-14 DEAL::Writing grid to <grid-4.eps> DEAL::Writing solution to <sol-4.gnuplot> DEAL::Cycle 5 DEAL::Number of active cells: 1561 DEAL::Number of degrees of freedom: 6244 DEAL:Richardson::Starting value 0.131369 DEAL:Richardson::Convergence step 81 value 2.39812e-13 DEAL::Writing grid to <grid-5.eps> DEAL::Writing solution to <sol-5.gnuplot>
We show the solutions on the initial mesh, the mesh after two and after five adaptive refinement steps.
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:
And finally we show a plot of a 3d computation.
(If you are looking at a locally installed deal.II version, then the program can be found at /home /hazelsct /repositories /deal.ii /examples /step-12 /step-12.cc . Otherwise, this is only the path on some remote server.)
/ * Subversion Id: @ref step_12 "step-12".cc 21222 2010-06-17 15:51:46Z kanschat * / / * Author: Guido Kanschat, Texas A&M University, 2009 * / / * Subversion Id: @ref step_12 "step-12".cc 21222 2010-06-17 15:51:46Z kanschat * / / * * / / * Copyright (C) 2010 by the deal.II authors * / / * * / / * This file is subject to QPL and may not be distributed * / / * without copyright and license information. Please refer * / / * to the file deal.II/doc/license.html for the text and * / / * further information on this license. * / #include <base/quadrature_lib.h> #include <base/function.h> #include <lac/vector.h> #include <lac/compressed_sparsity_pattern.h> #include <lac/sparse_matrix.h> #include <grid/tria.h> #include <grid/grid_generator.h> #include <grid/grid_out.h> #include <grid/grid_refinement.h> #include <grid/tria_accessor.h> #include <grid/tria_iterator.h> #include <fe/fe_values.h> #include <dofs/dof_handler.h> #include <dofs/dof_accessor.h> #include <dofs/dof_tools.h> #include <numerics/data_out.h> #include <fe/mapping_q1.h> #include <fe/fe_dgq.h> #include <lac/solver_richardson.h> #include <lac/precondition_block.h> #include <numerics/derivative_approximation.h> #include <numerics/mesh_worker.h> #include <numerics/mesh_worker_info.h> #include <numerics/mesh_worker_assembler.h> #include <numerics/mesh_worker_loop.h> #include <iostream> #include <fstream> using namespace dealii;
template <int dim> class BoundaryValues: public Function<dim> { public: BoundaryValues () {}; virtual void value_list (const std::vector<Point<dim> > &points, std::vector<double> &values, const unsigned int component=0) const; }; template <int dim> void BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points, std::vector<double> &values, const unsigned int) const { 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.; } }
template <int dim> class Step12 { public: Step12 (); 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; FE_DGQ<dim> fe; DoFHandler<dim> dof_handler; SparsityPattern sparsity_pattern; SparseMatrix<double> system_matrix; Vector<double> solution; Vector<double> right_hand_side; typedef MeshWorker::DoFInfo<dim> DoFInfo; typedef MeshWorker::IntegrationInfo<dim> CellInfo; 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> Step12<dim>::Step12 () : fe (1), dof_handler (triangulation) {} template <int dim> void Step12<dim>::setup_system () { dof_handler.distribute_dofs (fe); CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); DoFTools::make_flux_sparsity_pattern (dof_handler, c_sparsity); sparsity_pattern.copy_from(c_sparsity); system_matrix.reinit (sparsity_pattern); solution.reinit (dof_handler.n_dofs()); right_hand_side.reinit (dof_handler.n_dofs()); }
template <int dim> void Step12<dim>::assemble_system () { MeshWorker::IntegrationInfoBox<dim> info_box; 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); 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); info_box.initialize(fe, mapping); MeshWorker::DoFInfo<dim> dof_info(dof_handler); MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double> > assembler; assembler.initialize(system_matrix, right_hand_side); MeshWorker::integration_loop<dim, dim> (dof_handler.begin_active(), dof_handler.end(), dof_info, info_box, &Step12<dim>::integrate_cell_term, &Step12<dim>::integrate_boundary_term, &Step12<dim>::integrate_face_term, assembler, true); }
template <int dim> void Step12<dim>::integrate_cell_term (DoFInfo& dinfo, CellInfo& info) { const FEValuesBase<dim>& fe_v = info.fe_values(); FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix; const std::vector<double> &JxW = fe_v.get_JxW_values (); for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point) { Point<dim> beta; beta(0) = -fe_v.quadrature_point(point)(1); beta(1) = fe_v.quadrature_point(point)(0); beta /= beta.norm(); for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j) local_matrix(i,j) -= beta*fe_v.shape_grad(i,point)* fe_v.shape_value(j,point) * JxW[point]; } } template <int dim> void Step12<dim>::integrate_boundary_term (DoFInfo& dinfo, CellInfo& info) { const FEValuesBase<dim>& fe_v = 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_v.get_JxW_values (); const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors (); std::vector<double> g(fe_v.n_quadrature_points); static BoundaryValues<dim> boundary_function; boundary_function.value_list (fe_v.get_quadrature_points(), g); for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point) { Point<dim> beta; beta(0) = -fe_v.quadrature_point(point)(1); beta(1) = fe_v.quadrature_point(point)(0); beta /= beta.norm(); const double beta_n=beta * normals[point]; if (beta_n>0) for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j) local_matrix(i,j) += beta_n * fe_v.shape_value(j,point) * fe_v.shape_value(i,point) * JxW[point]; else for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) local_vector(i) -= beta_n * g[point] * fe_v.shape_value(i,point) * JxW[point]; } } template <int dim> void Step12<dim>::integrate_face_term (DoFInfo& dinfo1, DoFInfo& dinfo2, CellInfo& info1, CellInfo& info2) { const FEValuesBase<dim>& fe_v = info1.fe_values(); const FEValuesBase<dim>& fe_v_neighbor = info2.fe_values(); 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; const std::vector<double> &JxW = fe_v.get_JxW_values (); const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors (); for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point) { Point<dim> beta; beta(0) = -fe_v.quadrature_point(point)(1); beta(1) = fe_v.quadrature_point(point)(0); beta /= beta.norm(); const double beta_n=beta * normals[point]; if (beta_n>0) { for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j) u1_v1_matrix(i,j) += beta_n * fe_v.shape_value(j,point) * fe_v.shape_value(i,point) * JxW[point]; for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k) for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j) u1_v2_matrix(k,j) -= beta_n * fe_v.shape_value(j,point) * fe_v_neighbor.shape_value(k,point) * JxW[point]; } else { for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l) u2_v1_matrix(i,l) += beta_n * fe_v_neighbor.shape_value(l,point) * fe_v.shape_value(i,point) * JxW[point]; for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k) for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l) u2_v2_matrix(k,l) -= beta_n * fe_v_neighbor.shape_value(l,point) * fe_v_neighbor.shape_value(k,point) * JxW[point]; } } }
template <int dim> void Step12<dim>::solve (Vector<double> &solution) { SolverControl solver_control (1000, 1e-12); SolverRichardson<> solver (solver_control); PreconditionBlockSSOR<SparseMatrix<double> > preconditioner; preconditioner.initialize(system_matrix, fe.dofs_per_cell); solver.solve (system_matrix, solution, right_hand_side, preconditioner); } template <int dim> void Step12<dim>::refine_grid () { Vector<float> gradient_indicator (triangulation.n_active_cells()); DerivativeApproximation::approximate_gradient (mapping, dof_handler, solution, gradient_indicator); typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) gradient_indicator(cell_no)*=std::pow(cell->diameter(), 1+1.0*dim/2); GridRefinement::refine_and_coarsen_fixed_number (triangulation, gradient_indicator, 0.3, 0.1); triangulation.execute_coarsening_and_refinement (); } template <int dim> void Step12<dim>::output_results (const unsigned int cycle) const { std::string filename = "grid-"; filename += ('0' + cycle); Assert (cycle < 10, ExcInternalError()); filename += ".eps"; deallog << "Writing grid to <" << filename << ">" << std::endl; std::ofstream eps_output (filename.c_str()); GridOut grid_out; grid_out.write_eps (triangulation, eps_output); filename = "sol-"; filename += ('0' + cycle); Assert (cycle < 10, ExcInternalError()); filename += ".gnuplot"; deallog << "Writing solution to <" << filename << ">" << std::endl; std::ofstream gnuplot_output (filename.c_str()); 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); } template <int dim> void Step12<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); } } int main () { try { Step12<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; }