// The next new include file declares a base class <code>TensorFunction</code>
// not unlike the <code>Function</code> class, but with the difference that
-// the return value is tensor-valued rather than scalar of vector-valued.
+// TensorFunction::value returns a Tensor instead of a scalar.
#include <deal.II/base/tensor_function.h>
#include <deal.II/numerics/error_estimator.h>
// @sect3{AdvectionProblem class declaration}
- // Following we declare the main class of this program. It is very much
- // like the main classes of previous examples, so we again only comment on
- // the differences.
+ // Here comes the main class of this program. It is very much like the main
+ // classes of previous examples, so we again only comment on the
+ // differences.
template <int dim>
class AdvectionProblem
{
void copy_local_to_global(const AssemblyCopyData ©_data);
- // The following functions again are as in previous examples, as are the
- // subsequent variables.
+ // The following functions again are the same as they were in previous
+ // examples, as are the subsequent variables:
void solve();
void refine_grid();
void output_results(const unsigned int cycle) const;
// <code>Function</code> base class, as we have done for boundary values and
// coefficients in previous examples, but there is another possibility in
// the library, namely a base class that describes tensor valued
- // functions. In contrast to the usual <code>Function</code> objects, we
- // provide the compiler with knowledge on the size of the objects of the
- // return type. This enables the compiler to generate efficient code, which
- // is not so simple for usual vector-valued functions where memory has to be
- // allocated on the heap (thus, the <code>Function::vector_value</code>
- // function has to be given the address of an object into which the result
- // is to be written, in order to avoid copying and memory allocation and
- // deallocation on the heap). In addition to the known size, it is possible
- // not only to return vectors, but also tensors of higher rank; however,
- // this is not very often requested by applications, to be honest...
- //
- // The interface of the <code>TensorFunction</code> class is relatively
- // close to that of the <code>Function</code> class, so there is probably no
- // need to comment in detail the following declaration:
+ // functions. This is more convenient than overriding Function::value() with
+ // a method that knows about multiple function components: at the end of the
+ // day we need a Tensor, so we may as well just use a class that returns a
+ // Tensor.
template <int dim>
class AdvectionField : public TensorFunction<1, dim>
{
// respective macro <code>DeclExceptionN</code>.
//
// To learn how the preprocessor expands this macro into actual code,
- // please refer to the documentation of the exception classes in the base
- // library. Suffice it to say that by this macro call, the respective
- // exception class is declared, which also has error output functions
- // already implemented.
+ // please refer to the documentation of the exception classes. In brief,
+ // this macro call declares and defines a class
+ // <code>ExcDimensionMismatch</code> inheriting from ExceptionBase which
+ // implements all necessary error output functions.
};
-
-
// The following two functions implement the interface described above. The
// first simply implements the function as described in the introduction,
// while the second uses the same trick to avoid calling a virtual function
Point<dim> value;
value[0] = 2;
for (unsigned int i = 1; i < dim; ++i)
- value[i] = 1 + 0.8 * std::sin(8 * numbers::PI * p[0]);
+ value[i] = 1 + 0.8 * std::sin(8. * numbers::PI * p[0]);
return value;
}
// Besides the advection field, we need two functions describing the source
- // terms (<code>right hand side</code>) and the boundary values. First for
- // the right hand side, which follows the same pattern as in previous
- // examples. As described in the introduction, the source is a constant
- // function in the vicinity of a source point, which we denote by the
- // constant static variable <code>center_point</code>. We set the values of
- // this center using the same template tricks as we have shown in the step-7
- // example program. The rest is simple and has been shown previously.
+ // terms (<code>right hand side</code>) and the boundary values. As
+ // described in the introduction, the source is a constant function in the
+ // vicinity of a source point, which we denote by the constant static
+ // variable <code>center_point</code>. We set the values of this center
+ // using the same template tricks as we have shown in the step-7 example
+ // program. The rest is simple and has been shown previously.
template <int dim>
class RightHandSide : public Function<dim>
{
Assert(component == 0, ExcIndexRange(component, 0, 1));
const double diameter = 0.1;
return ((p - center_point).norm_square() < diameter * diameter ?
- .1 / std::pow(diameter, dim) :
- 0);
+ 0.1 / std::pow(diameter, dim) :
+ 0.0);
}
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
- const double sine_term =
- std::sin(16 * numbers::PI * std::sqrt(p.norm_square()));
- const double weight = std::exp(-5 * p.norm_square()) / std::exp(-5.);
+ const double sine_term = std::sin(16. * numbers::PI * p.norm());
+ const double weight = std::exp(-5. * p.norm_square()) / std::exp(-5.);
return sine_term * weight;
}
// Now, finally, here comes the class that will compute the difference
// approximation of the gradient on each cell and weighs that with a power
- // of the mesh size, as described in the introduction. This class is a
+ // of the mesh size, as described in the introduction. This class is a
// simple version of the <code>DerivativeApproximation</code> class in the
// library, that uses similar techniques to obtain finite difference
// approximations of the gradient of a finite element field, or of higher
// to the function rather than to the class itself. This way, you don't have
// to specify the template parameter yourself as in most other cases, but
// the compiler can figure its value out itself from the dimension of the
- // DoF handler object that one passes as first argument.
+ // DoFHandler object that one passes as first argument.
//
// Before jumping into the fray with the implementation, let us also comment
// on the parallelization strategy. We have already introduced the necessary
// framework for using the WorkStream concept in the declaration of the main
// class of this program above. We will use it again here. In the current
- // context, this means that we have to define (i) classes for scratch and
- // copy objects, (ii) a function that does the local computation on one
- // cell, and (iii) a function that copies the local result into a global
- // object. Given this general framework, we will, however, deviate from it a
+ // context, this means that we have to define
+ // <ol>
+ // <li>classes for scratch and copy objects,</li>
+ // <li>a function that does the local computation on one cell, and</li>
+ // <li>a function that copies the local result into a global object.</li>
+ // </ol>
+ // Given this general framework, we will, however, deviate from it a
// bit. In particular, WorkStream was generally invented for cases where
// each local computation on a cell <i>adds</i> to a global object -- for
// example, when assembling linear systems where we add local contributions
// into a global matrix and right hand side. WorkStream is designed to handle
// the potential conflict of multiple threads trying to do this addition at
// the same time, and consequently has to provide for some way to ensure that
- // only thread gets to do this at a time. Here, however, the situation is
+ // only one thread gets to do this at a time. Here, however, the situation is
// slightly different: we compute contributions from every cell
// individually, but then all we need to do is put them into an element of
// an output vector that is unique to each cell. Consequently, there is no
// risk that the write operations from two cells might conflict, and the
// elaborate machinery of WorkStream to avoid conflicting writes is not
// necessary. Consequently, what we will do is this: We still need a scratch
- // object that holds, for example, the FEValues object. However, we only
+ // object that holds, for example, the FEValues object. However, we only
// create a fake, empty copy data structure. Likewise, we do need the
// function that computes local contributions, but since it can already put
// the result into its final location, we do not need a copy-local-to-global
// contain FEValues and FEFaceValues objects, and so we will need to
// have constructors and copy constructors that allow us to create
// them. In initializing them, note first that we use bilinear
- // elements, soGauss formulae with two points in each space
- // direction are sufficient. For the cell terms we need the values
+ // elements, so Gauss formulae with two points in each space
+ // direction are sufficient. For the cell terms we need the values
// and gradients of the shape functions, the quadrature points in
// order to determine the source density and the advection field at
// a given point, and the weights of the quadrature points times the
// Now, this is the function that does the actual work. It is not very
// different from the <code>assemble_system</code> functions of previous
// example programs, so we will again only comment on the differences. The
- // mathematical stuff follows closely what we have said in the introduction.
+ // mathematical stuff closely follows what we have said in the introduction.
//
// There are a number of points worth mentioning here, though. The
// first one is that we have moved the FEValues and FEFaceValues
// the direction of flow at this point; we obtain this information
// using the FEFaceValues object and only decide within the main loop
// whether a quadrature point is on the inflow boundary.
- for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
- ++face)
- if (cell->face(face)->at_boundary())
+ for (unsigned int face_n = 0; face_n < GeometryInfo<dim>::faces_per_cell;
+ ++face_n)
+ if (cell->face(face_n)->at_boundary())
{
// Ok, this face of the present cell is on the boundary of the
// domain. Just as for the usual FEValues object which we have
// used in previous examples and also above, we have to
// reinitialize the FEFaceValues object for the present face:
- scratch_data.fe_face_values.reinit(cell, face);
+ scratch_data.fe_face_values.reinit(cell, face_n);
// For the quadrature points at hand, we ask for the values of
// the inflow function and for the direction of flow:
scratch_data.fe_face_values.get_quadrature_points(),
face_advection_directions);
- // Now loop over all quadrature points and see whether it is on
- // the inflow or outflow part of the boundary. This is
- // determined by a test whether the advection direction points
- // inwards or outwards of the domain (note that the normal
- // vector points outwards of the cell, and since the cell is at
- // the boundary, the normal vector points outward of the domain,
+ // Now loop over all quadrature points and see whether this face is on
+ // the inflow or outflow part of the boundary. The normal
+ // vector points out of the cell: since the face is at
+ // the boundary, the normal vector points out of the domain,
// so if the advection direction points into the domain, its
- // scalar product with the normal vector must be negative):
+ // scalar product with the normal vector must be negative (to see why
+ // this is true, consider the scalar product definition that uses a
+ // cosine):
for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)
if (scratch_data.fe_face_values.normal_vector(q_point) *
face_advection_directions[q_point] <
- 0)
- // If the is part of the inflow boundary, then compute the
+ 0.)
+ // If the face is part of the inflow boundary, then compute the
// contributions of this face to the global matrix and right
// hand side, using the values obtained from the
// FEFaceValues object and the formulae discussed in the
}
}
-
- // Now go on by transferring the local contributions to the system of
- // equations into the global objects. The first step was to obtain the
- // global indices of the degrees of freedom on this cell.
+ // The final piece of information the copy routine needs is the global
+ // indices of the degrees of freedom on this cell, so we end by writing
+ // them to the local array:
cell->get_dof_indices(copy_data.local_dof_indices);
}
// The second function we needed to write was the one that copies
- // the local contributions the previous function has computed and
- // put into the copy data object, into the global matrix and right
+ // the local contributions the previous function computed (and
+ // put into the AssemblyCopyData object) into the global matrix and right
// hand side vector objects. This is essentially what we always had
// as the last block of code when assembling something on every
// cell. The following should therefore be pretty obvious:
}
}
-
-
- // Following is the function that solves the linear system of equations. As
- // the system is no more symmetric positive definite as in all the previous
- // examples, we can't use the Conjugate Gradients method anymore. Rather, we
- // use a solver that is tailored to nonsymmetric systems like the one at
- // hand, the BiCGStab method. As preconditioner, we use the Jacobi method.
+ // Here comes the linear solver routine. As the system is no longer
+ // symmetric positive definite as in all the previous examples, we cannot
+ // use the Conjugate Gradient method anymore. Rather, we use a solver that
+ // is tailored to nonsymmetric systems like the one at hand, the BiCGStab
+ // method. As preconditioner, we use the Jacobi method.
template <int dim>
void AdvectionProblem<dim>::solve()
{
hanging_node_constraints.distribute(solution);
}
-
// The following function refines the grid according to the quantity
// described in the introduction. The respective computations are made in
// the class <code>GradientEstimation</code>. The only difference to
triangulation.execute_coarsening_and_refinement();
}
-
-
// Writing output to disk is done in the same way as in the previous
// examples. Indeed, the function is identical to the one in step-6.
template <int dim>
}
- // ... as is the main loop (setup -- solve -- refine)
+ // ... as is the main loop (setup -- solve -- refine), aside from the number
+ // of cycles and the initial grid:
template <int dim>
void AdvectionProblem<dim>::run()
{
{}
- // Next for the implementation of the <code>GradientEstimation</code>
+ // Next comes the implementation of the <code>GradientEstimation</code>
// class. The first function does not much except for delegating work to the
// other function, but there is a bit of setup at the top.
//
// Before starting with the work, we check that the vector into
// which the results are written has the right size. Programming
// mistakes in which one forgets to size arguments correctly at the
- // calling site are quite common. Because the resulting damage from
+ // calling site are quite common. Because the resulting damage from
// not catching such errors is often subtle (e.g., corruption of
// data somewhere in memory, or non-reproducible results), it is
// well worth the effort to check for such things.
}
- // Following now the function that actually computes the finite difference
- // approximation to the gradient. The general outline of the function is to
- // first compute the list of active neighbors of the present cell and then
- // compute the quantities described in the introduction for each of the
- // neighbors. The reason for this order is that it is not a one-liner to
- // find a given neighbor with locally refined meshes. In principle, an
+ // Here comes the function that estimates the local error by computing the
+ // finite difference approximation of the gradient. The function first
+ // computes the list of active neighbors of the present cell and then
+ // computes the quantities described in the introduction for each of
+ // the neighbors. The reason for this order is that it is not a one-liner
+ // to find a given neighbor with locally refined meshes. In principle, an
// optimized implementation would find neighbors and the quantities
// depending on them in one step, rather than first building a list of
// neighbors and in a second step their contributions but we will gladly
// passed to WorkStream::run works on "scratch" objects that keep all
// temporary objects. This way, we do not need to create and initialize
// objects that are expensive to initialize within the function that does
- // the work, every time it is called for a given cell. Such an argument is
+ // the work every time it is called for a given cell. Such an argument is
// passed as the second argument. The third argument would be a "copy-data"
// object (see @ref threads for more information) but we do not actually use
- // any of these here. Because WorkStream::run() insists on passing three
+ // any of these here. Since WorkStream::run() insists on passing three
// arguments, we declare this function with three arguments, but simply
// ignore the last one.
//
- // (This is unsatisfactory from an esthetic perspective. It can be avoided,
+ // (This is unsatisfactory from an aesthetic perspective. It can be avoided,
// at the cost of some other trickery. If you allow, let us here show
// how. First, assume that we had declared this function to only take two
// arguments by omitting the unused last one. Now, WorkStream::run still
// way to "forget" the third argument in the call. Simply passing
// WorkStream::run the pointer to the function as we do above will not do
// this -- the compiler will complain that a function declared to have two
- // arguments is called with three arguments. However, we can do this by
+ // arguments is called with three arguments. However, we can do this by
// passing the following as the third argument when calling WorkStream::run()
// above:
// @code
// <code>Y</code> tensor:
scratch_data.fe_midpoint_value.reinit(cell);
- // Then allocate the vector that will be the sum over the y-vectors
- // times the approximate directional derivative:
- Tensor<1, dim> projected_gradient;
-
-
- // Now before going on first compute a list of all active neighbors of
- // the present cell. We do so by first looping over all faces and see
+ // Now, before we go on, we first compute a list of all active neighbors
+ // of the present cell. We do so by first looping over all faces and see
// whether the neighbor there is active, which would be the case if it
// is on the same level as the present cell or one level coarser (note
// that a neighbor can only be once coarser than the present cell, as
// have to clear the array storing the iterators to the active
// neighbors, of course.
active_neighbors.clear();
- for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
- ++face_no)
- if (!cell->at_boundary(face_no))
+ for (unsigned int face_n = 0; face_n < GeometryInfo<dim>::faces_per_cell;
+ ++face_n)
+ if (!cell->at_boundary(face_n))
{
// First define an abbreviation for the iterator to the face and
// the neighbor
const typename DoFHandler<dim>::face_iterator face =
- cell->face(face_no);
+ cell->face(face_n);
const typename DoFHandler<dim>::cell_iterator neighbor =
- cell->neighbor(face_no);
+ cell->neighbor(face_n);
// Then check whether the neighbor is active. If it is, then it
// is on the same level or one level coarser (if we are not in
typename DoFHandler<dim>::cell_iterator neighbor_child =
neighbor;
while (neighbor_child->has_children())
- neighbor_child =
- neighbor_child->child(face_no == 0 ? 1 : 0);
+ neighbor_child = neighbor_child->child(face_n == 0 ? 1 : 0);
// As this used some non-trivial geometrical intuition,
// we might want to check whether we did it right,
- // i.e. check whether the neighbor of the cell we found
+ // i.e., check whether the neighbor of the cell we found
// is indeed the cell we are presently working
// on. Checks like this are often useful and have
// frequently uncovered errors both in algorithms like
// an error that is irrecoverable and probably qualifies
// as an internal error. We therefore use a predefined
// exception class to throw here.
- Assert(neighbor_child->neighbor(face_no == 0 ? 1 : 0) == cell,
+ Assert(neighbor_child->neighbor(face_n == 0 ? 1 : 0) == cell,
ExcInternalError());
// If the check succeeded, we push the active neighbor
}
else
// If we are not in 1d, we collect all neighbor children
- // `behind' the subfaces of the current face
- for (unsigned int subface_no = 0;
- subface_no < face->n_children();
- ++subface_no)
+ // `behind' the subfaces of the current face and move on:
+ for (unsigned int subface_n = 0; subface_n < face->n_children();
+ ++subface_n)
active_neighbors.push_back(
- cell->neighbor_child_on_subface(face_no, subface_no));
+ cell->neighbor_child_on_subface(face_n, subface_n));
}
}
scratch_data.fe_midpoint_value.get_function_values(scratch_data.solution,
this_midpoint_value);
-
// Now loop over all active neighbors and collect the data we
// need. Allocate a vector just like <code>this_midpoint_value</code>
// which we will use to store the value of the solution in the
// iteration of this inner loop (memory allocation is a rather
// expensive operation):
std::vector<double> neighbor_midpoint_value(1);
+ Tensor<1, dim> projected_gradient;
typename std::vector<typename DoFHandler<dim>::active_cell_iterator>::
const_iterator neighbor_ptr = active_neighbors.begin();
for (; neighbor_ptr != active_neighbors.end(); ++neighbor_ptr)
*neighbor_ptr;
// Then get the center of the neighbor cell and the value of the
- // finite element function thereon. Note that for this information
- // we have to reinitialize the <code>FEValues</code> object for
- // the neighbor cell.
+ // finite element function at that point. Note that for this
+ // information we have to reinitialize the <code>FEValues</code>
+ // object for the neighbor cell.
scratch_data.fe_midpoint_value.reinit(neighbor);
const Point<dim> neighbor_center =
scratch_data.fe_midpoint_value.quadrature_point(0);
// span the whole space, otherwise we would not have all components of
// the gradient. This is indicated by the invertibility of the matrix.
//
- // If the matrix should not be invertible, this means that the present
+ // If the matrix is not invertible, then the present
// cell had an insufficient number of active neighbors. In contrast to
- // all previous cases, where we raised exceptions, this is, however,
+ // all previous cases (where we raised exceptions) this is, however,
// not a programming error: it is a runtime error that can happen in
// optimized mode even if it ran well in debug mode, so it is
// reasonable to try to catch this error also in optimized mode. For
// this case, there is the <code>AssertThrow</code> macro: it checks
// the condition like the <code>Assert</code> macro, but not only in
// debug mode; it then outputs an error message, but instead of
- // terminating the program as in the case of the <code>Assert</code>
+ // aborting the program as in the case of the <code>Assert</code>
// macro, the exception is thrown using the <code>throw</code> command
// of C++. This way, one has the possibility to catch this error and
// take reasonable counter actions. One such measure would be to
// least once.
AssertThrow(determinant(Y) != 0, ExcInsufficientDirections());
- // If, on the other hand the matrix is invertible, then invert it,
- // multiply the other quantity with it and compute the estimated error
- // using this quantity and the right powers of the mesh width:
+ // If, on the other hand, the matrix is invertible, then invert it,
+ // multiply the other quantity with it, and compute the estimated error
+ // using this quantity and the correct powers of the mesh width:
const Tensor<2, dim> Y_inverse = invert(Y);
- Tensor<1, dim> gradient = Y_inverse * projected_gradient;
+ const Tensor<1, dim> gradient = Y_inverse * projected_gradient;
// The last part of this function is the one where we write into
// the element of the output vector what we have just
// at the correct element inside this vector -- but we can ask the
// cell we're on the how-manyth active cell it is for this:
scratch_data.error_per_cell(cell->active_cell_index()) =
- (std::pow(cell->diameter(), 1 + 1.0 * dim / 2) *
- std::sqrt(gradient.norm_square()));
+ (std::pow(cell->diameter(), 1 + 1.0 * dim / 2) * gradient.norm());
}
} // namespace Step9
// @sect3{Main function}
-// The <code>main</code> function is similar to the previous examples. The main
-// difference is that we use MultithreadInfo to set the maximum
-// number of threads (see @ref threads "Parallel computing with multiple
-// processors accessing shared memory" documentation module for more
-// explanation). The number of threads used is the minimum of the environment
-// variable DEAL_II_NUM_THREADS and the parameter of
+// The <code>main</code> function is similar to the previous examples. The
+// primary difference is that we use MultithreadInfo to set the maximum
+// number of threads (see the documentation module @ref threads
+// "Parallel computing with multiple processors accessing shared memory"
+// for more information). The number of threads used is the minimum of the
+// environment variable DEAL_II_NUM_THREADS and the parameter of
// <code>set_thread_limit</code>. If no value is given to
// <code>set_thread_limit</code>, the default value from the Intel Threading
// Building Blocks (TBB) library is used. If the call to
// <code>set_thread_limit</code> is omitted, the number of threads will be
-// chosen by TBB indepently of DEAL_II_NUM_THREADS.
+// chosen by TBB independently of DEAL_II_NUM_THREADS.
int main()
{
+ using namespace dealii;
try
{
- dealii::MultithreadInfo::set_thread_limit();
+ MultithreadInfo::set_thread_limit();
Step9::AdvectionProblem<2> advection_problem_2d;
advection_problem_2d.run();
projects / dealii.git / commitdiff