From 474b7066202ca067ef3785f83e6af721e90dda20 Mon Sep 17 00:00:00 2001 From: bangerth Date: Sat, 23 Nov 2013 01:57:28 +0000 Subject: [PATCH] Document the most recent WorkStream changes to step-9. Update output. Use VTK format instead of the outdated GMV format. git-svn-id: https://svn.dealii.org/trunk@31769 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-9/doc/results.dox | 18 +-- deal.II/examples/step-9/step-9.cc | 167 ++++++++++++++---------- 2 files changed, 106 insertions(+), 79 deletions(-) diff --git a/deal.II/examples/step-9/doc/results.dox b/deal.II/examples/step-9/doc/results.dox index 13ff02506b..d0e541a513 100644 --- a/deal.II/examples/step-9/doc/results.dox +++ b/deal.II/examples/step-9/doc/results.dox @@ -12,17 +12,17 @@ Cycle 1: Number of active cells: 643 Number of degrees of freedom: 793 Cycle 2: - Number of active cells: 1663 - Number of degrees of freedom: 1946 + Number of active cells: 1669 + Number of degrees of freedom: 1950 Cycle 3: - Number of active cells: 4219 - Number of degrees of freedom: 4905 + Number of active cells: 4231 + Number of degrees of freedom: 4923 Cycle 4: - Number of active cells: 10708 - Number of degrees of freedom: 12128 + Number of active cells: 10753 + Number of degrees of freedom: 12175 Cycle 5: - Number of active cells: 26908 - Number of degrees of freedom: 29698 + Number of active cells: 27004 + Number of degrees of freedom: 29810 @endcode As can be seen, quite a number of cells is used on the finest level to resolve the features of the solution. The final grid showing this is @@ -42,7 +42,7 @@ solution itself: Note that the solution is created by that part that is transported -along the wiggled advection field from the left and lower boundaries +along the wiggly advection field from the left and lower boundaries to the top right, and the part that is created by the source in the lower left corner, and the results of which are also transported along. The grid shown above is well-adapted to resolve these diff --git a/deal.II/examples/step-9/step-9.cc b/deal.II/examples/step-9/step-9.cc index c10edf29eb..c1be309814 100644 --- a/deal.II/examples/step-9/step-9.cc +++ b/deal.II/examples/step-9/step-9.cc @@ -961,18 +961,21 @@ namespace Step9 data_out.add_data_vector (solution, "solution"); data_out.build_patches (); - std::ofstream output ("final-solution.gmv"); - data_out.write_gmv (output); + std::ofstream output ("final-solution.vtk"); + data_out.write_vtk (output); } // @sect3{GradientEstimation class implementation} - // ScratchData used by estimate_cell + // Now for the implementation of the GradientEstimation class. + // Let us start by defining constructors for the + // EstimateScratchData class used by the + // estimate_cell() function: template - GradientEstimation::EstimateScratchData - ::EstimateScratchData (const FiniteElement &fe, + GradientEstimation::EstimateScratchData:: + EstimateScratchData (const FiniteElement &fe, const Vector &solution) : fe_midpoint_value(fe, @@ -982,11 +985,9 @@ namespace Step9 {} - - // ScratchData used by estimate_cell template - GradientEstimation::EstimateScratchData - ::EstimateScratchData(const EstimateScratchData &scratch_data) + GradientEstimation::EstimateScratchData:: + EstimateScratchData(const EstimateScratchData &scratch_data) : fe_midpoint_value(scratch_data.fe_midpoint_value.get_fe(), scratch_data.fe_midpoint_value.get_quadrature(), @@ -995,88 +996,108 @@ namespace Step9 {} - // Now for the implementation of the GradientEstimation + // Next for the implementation of the GradientEstimation // class. The first function does not much except for delegating work to the - // other function: + // 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. It is a common error that such + // parameters have the wrong size, but the resulting damage by not + // catching these errors are very subtle as they are usually corruption of + // data somewhere in memory. Often, the problems emerging from this are + // not reproducible, and it is well worth the effort to + // check for such things. + // + // The second piece is to set up the iterator that goes in lockstep over the + // cells of the domain and the corresponding elements of the output vector + // (see above where we introduced the SynchronousIterators + // class). We can abbreviate the process slightly by introducing a + // typedef that denotes a pair of iterators. This being set up, + // we can hand the whole thing off to WorkStream::run, keeping in mind that + // we do not need a copy-local-to-global function here but can get away by + // simply using a default-constructed function object (the equivalent to a + // NULL function pointer). template void GradientEstimation::estimate (const DoFHandler &dof_handler, const Vector &solution, Vector &error_per_cell) { - // Before starting with the work, we check that the vector into which the - // results are written, has the right size. It is a common error that such - // parameters have the wrong size, but the resulting damage by not - // catching these errors are very subtle as they are usually corruption of - // data somewhere in memory. Often, the problems emerging from this are - // not reproducible, and we found that it is well worth the effort to - // check for such things. Assert (error_per_cell.size() == dof_handler.get_tria().n_active_cells(), ExcInvalidVectorLength (error_per_cell.size(), dof_handler.get_tria().n_active_cells())); - // In the same way as before, we use a Threads::ThreadGroup - // object to collect the descriptor objects of different threads. Note - // that as the function called is not a member function, but rather a - // static function, we need not (and can not) pass a this - // pointer to the new_thread function in this case. - // - // Taking pointers to templated functions seems to be notoriously - // difficult for many compilers (since there are several functions with - // the same name -- just as with overloaded functions). It therefore - // happens quite frequently that we can't directly insert taking the - // address of a function in the call to encapsulate for one - // or the other compiler, but have to take a temporary variable for that - // purpose. Here, in this case, Compaq's cxx compiler choked - // on the code so we use this workaround with the function pointer: - void (*estimate_cell_ptr) (const SynchronousIterators::active_cell_iterator,Vector::iterator> > &cell, - EstimateScratchData &scratch_data, - const EstimateCopyData ©_data) - = &GradientEstimation::template estimate_cell; - typedef std_cxx1x::tuple::active_cell_iterator,Vector::iterator> - Iterators; - SynchronousIterators begin_sync_it(Iterators(dof_handler.begin_active(), - error_per_cell.begin())); - SynchronousIterators end_sync_it(Iterators(dof_handler.end(),error_per_cell.end())); - - WorkStream::run(begin_sync_it,end_sync_it, - estimate_cell_ptr, - std_cxx1x::function (), - EstimateScratchData (dof_handler.get_fe(),solution), - EstimateCopyData ()); - - // Note that if the value of the variable - // multithread_info.n_threads() was one, or if the - // library was not configured to use threads, then the sequence of - // commands above reduced to a complicated way to simply call the - // estimate_interval function with the whole range of cells - // to work on. However, using the way above, we are able to write the - // program such that it makes no difference whether we presently work with - // multiple threads or in single-threaded mode, thus eliminating the need - // to write code included in conditional preprocessor sections. + IteratorTuple; + + SynchronousIterators + begin_sync_it (IteratorTuple (dof_handler.begin_active(), + error_per_cell.begin())), + end_sync_it (IteratorTuple (dof_handler.end(), + error_per_cell.end())); + + WorkStream::run (begin_sync_it, + end_sync_it, + &GradientEstimation::template estimate_cell, + std_cxx1x::function (), + EstimateScratchData (dof_handler.get_fe(), + solution), + EstimateCopyData ()); } // Following now the function that actually computes the finite difference // approximation to the gradient. The general outline of the function is to - // loop over all the cells in the range of iterators designated by the third - // argument, and on each cell 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 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. + // 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 + // 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 + // leave this as an exercise. As discussed before, the worker function + // 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 + // 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 + // 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, + // 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 + // wants to call this function with three arguments, so we need to find a + // 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 + // passing the following as the third argument when calling WorkStream::run + // above: + // @code + // std_cxx1x::function &, + // EstimateScratchData &, + // EstimateCopyData &)> + // (std_cxx1x::bind (&GradientEstimation::template estimate_cell, + // std_cxx1x::_1, + // std_cxx1x::_2)) + // @endcode + // This creates a function object taking three arguments, but when it calls + // the underlying function object, it simply only uses the first and second + // argument -- we simply "forget" to use the third argument :-) + // In the end, this isn't completely obvious either, and so we didn't implement + // it, but hey -- it can be done!) // // Now for the details: template void - GradientEstimation::estimate_cell (const SynchronousIterators::active_cell_iterator,Vector::iterator> > &cell, - EstimateScratchData &scratch_data, - const EstimateCopyData ©_data) + GradientEstimation::estimate_cell (const SynchronousIterators::active_cell_iterator, + Vector::iterator> > &cell, + EstimateScratchData &scratch_data, + const EstimateCopyData &) { // We need space for the tensor Y, which is the sum of // outer products of the y-vectors. @@ -1284,6 +1305,12 @@ namespace Step9 Point gradient; contract (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 computed. As above, we need to get + // at the second element of the pair of iterators, which requires + // slightly awkward syntax but is not otherwise particularly + // difficult: *(std_cxx1x::get<1>(cell.iterators)) = (std::pow(std_cxx1x::get<0>(cell.iterators)->diameter(), 1+1.0*dim/2) * std::sqrt(gradient.square())); -- 2.39.5