/* $Id$ */
-/* Authors: Ralf Hartmann, University of Heidelberg, 2001 */
-/* hp-Version, Oliver Kayser-Herold, TU-Braunschweig 2005 */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2001, 2002 */
/* $Id$ */
/* Version: $Name$ */
/* */
-/* Copyright (C) 2001-2006 by the deal.II authors */
+/* Copyright (C) 2001, 2002, 2003, 2004, 2006 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. */
- // The first few files have already
- // been covered in previous examples
- // and will thus not be further
- // commented on.
+
+ // As in all programs, we start with
+ // a list of include files from the
+ // library, and as usual they are in
+ // the standard order which is
+ // <code>base</code> -- <code>lac</code> -- <code>grid</code> --
+ // <code>dofs</code> -- <code>fe</code> -- <code>numerics</code>
+ // (as each of these categories
+ // roughly builds upon previous
+ // ones), then C++ standard headers:
#include <base/quadrature_lib.h>
#include <base/function.h>
+#include <base/logstream.h>
+#include <base/table_handler.h>
+#include <base/thread_management.h>
#include <lac/vector.h>
+#include <lac/full_matrix.h>
#include <lac/sparse_matrix.h>
-#include <lac/vector_memory.h>
-#include <lac/solver_gmres.h>
+#include <lac/solver_cg.h>
#include <lac/precondition.h>
-#include <lac/sparse_ilu.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>
- // Instead of the usual DoFHandler
- // the hp::DoFHandler class has to be
- // used to gain access to the
- // hp-Functionality. The hpDoFHandler
- // provides essentially the same
- // interface as the standard DoFHandler.
-#include <dofs/hp_dof_handler.h>
- // Due to the implementation of the
- // hp-Method, the DoFAccessor classes
- // stay the same as for the DoFHandler,
- // and hence also require the same
- // include files.
+#include <grid/grid_refinement.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_constraints.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
-#include <dofs/dof_renumbering.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
#include <numerics/data_out.h>
+#include <numerics/error_estimator.h>
-#include <fe/mapping_q1.h>
-#include <fe/fe_dgq.h>
-
- // We are going to use gradients as
- // refinement indicator.
-#include <numerics/derivative_approximation.h>
- // Finally we do some time comparison
- // using the <code>Timer</code> class.
-#include <base/timer.h>
-
- // And this again is C++:
+ // Now for the C++ standard headers:
#include <iostream>
#include <fstream>
-
- // Now some additional utility classes will
- // be necessary. As the name "collection"
- // suggests, these are essentially
- // container classes which store the
- // quadrature and finite element objects,
- // for the different polynomial degrees.
-#include <fe/fe_collection.h>
-#include <fe/q_collection.h>
-
- // The first of the following
- // two files provides the hp::FEValues
- // class, which implements the same
- // functionality as the FEValues class
- // with the difference that it takes
- // "collection" objects instead of
- // single finite element or quadrature
- // objects. I.e. instead of the
- // usual Quadrature object a
- // QCollection object is needed.
-#include <fe/hp_fe_values.h>
-
- // A compressed sparsity pattern is
- // not an explicit prerequisite for the
- // use of the hp-functionality. But as
- // a standard sparsity pattern has to
- // based on a bandwidth estimate for the
- // highest polynomial degree, it would
- // be prohibitively bad. Therefore,
- // the recommended way is to explicitly
- // build a compressed sparsity pattern before
- // creating the matrices.
-#include <lac/compressed_sparsity_pattern.h>
-
+#include <list>
+#include <sstream>
// The last step is as in all
// previous programs:
using namespace dealii;
-
- // @sect3{Equation data}
+ // @sect3{Evaluation of the solution}
+
+ // As for the program itself, we
+ // first define classes that evaluate
+ // the solutions of a Laplace
+ // equation. In fact, they can
+ // evaluate every kind of solution,
+ // as long as it is described by a
+ // <code>DoFHandler</code> object, and a
+ // solution vector. We define them
+ // here first, even before the
+ // classes that actually generate the
+ // solution to be evaluated, since we
+ // need to declare an abstract base
+ // class that the solver classes can
+ // refer to.
//
- // First we define the classes
- // representing the equation-specific
- // functions. Both classes, <code>RHS</code>
- // and <code>BoundaryValues</code>, are
- // derived from the <code>Function</code>
- // class. Only the <code>value_list</code>
- // function are implemented because
- // only lists of function values are
- // computed rather than single
- // values.
-template <int dim>
-class RHS: public Function<dim>
+ // From an abstract point of view, we
+ // declare a pure base class
+ // that provides an evaluation
+ // operator <code>operator()</code> which will
+ // do the evaluation of the solution
+ // (whatever derived classes might
+ // consider an <code>evaluation</code>). Since
+ // this is the only real function of
+ // this base class (except for some
+ // bookkeeping machinery), one
+ // usually terms such a class that
+ // only has an <code>operator()</code> a
+ // <code>functor</code> in C++ terminology,
+ // since it is used just like a
+ // function object.
+ //
+ // Objects of this functor type will
+ // then later be passed to the solver
+ // object, which applies it to the
+ // solution just computed. The
+ // evaluation objects may then
+ // extract any quantity they like
+ // from the solution. The advantage
+ // of putting these evaluation
+ // functions into a separate
+ // hierarchy of classes is that by
+ // design they cannot use the
+ // internals of the solver object and
+ // are therefore independent of
+ // changes to the way the solver
+ // works. Furthermore, it is trivial
+ // to write another evaluation class
+ // without modifying the solver
+ // class, which speeds up programming
+ // (not being able to use internals
+ // of another class also means that
+ // you do not have to worry about
+ // them -- programming evaluators is
+ // usually a rather quickly done
+ // task), as well as compilation (if
+ // solver and evaluation classes are
+ // put into different files: the
+ // solver only needs to see the
+ // declaration of the abstract base
+ // class, and therefore does not need
+ // to be recompiled upon addition of
+ // a new evaluation class, or
+ // modification of an old one).
+ // On a related note, you can reuse
+ // the evaluation classes for other
+ // projects, solving different
+ // equations.
+ //
+ // In order to improve separation of
+ // code into different modules, we
+ // put the evaluation classes into a
+ // namespace of their own. This makes
+ // it easier to actually solve
+ // different equations in the same
+ // program, by assembling it from
+ // existing building blocks. The
+ // reason for this is that classes
+ // for similar purposes tend to have
+ // the same name, although they were
+ // developed in different
+ // contexts. In order to be able to
+ // use them together in one program,
+ // it is necessary that they are
+ // placed in different
+ // namespaces. This we do here:
+namespace Evaluation
{
- public:
- virtual void value_list (const std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component=0) const;
-};
+ // Now for the abstract base class
+ // of evaluation classes: its main
+ // purpose is to declare a pure
+ // virtual function <code>operator()</code>
+ // taking a <code>DoFHandler</code> object,
+ // and the solution vector. In
+ // order to be able to use pointers
+ // to this base class only, it also
+ // has to declare a virtual
+ // destructor, which however does
+ // nothing. Besides this, it only
+ // provides for a little bit of
+ // bookkeeping: since we usually
+ // want to evaluate solutions on
+ // subsequent refinement levels, we
+ // store the number of the present
+ // refinement cycle, and provide a
+ // function to change this number.
+ template <int dim>
+ class EvaluationBase
+ {
+ public:
+ virtual ~EvaluationBase ();
-template <int dim>
-class BoundaryValues: public Function<dim>
-{
- public:
- virtual void value_list (const std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component=0) const;
-};
+ void set_refinement_cycle (const unsigned int refinement_cycle);
+
+ virtual void operator () (const DoFHandler<dim> &dof_handler,
+ const Vector<double> &solution) const = 0;
+ protected:
+ unsigned int refinement_cycle;
+ };
+
+
+ // After the declaration has been
+ // discussed above, the
+ // implementation is rather
+ // straightforward:
+ template <int dim>
+ EvaluationBase<dim>::~EvaluationBase ()
+ {}
+
+
+ template <int dim>
+ void
+ EvaluationBase<dim>::set_refinement_cycle (const unsigned int step)
+ {
+ refinement_cycle = step;
+ }
- // The class <code>Beta</code> represents the
- // vector valued flow field of the
- // linear transport equation and is
- // not derived from the <code>Function</code>
- // class as we prefer to get function
- // values of type <code>Point</code> rather
- // than of type
- // <code>Vector@<double@></code>. This, because
- // there exist scalar products
- // between <code>Point</code> and <code>Point</code> as
- // well as between <code>Point</code> and
- // <code>Tensor</code>, simplifying terms like
- // $\beta\cdot n$ and
- // $\beta\cdot\nabla v$.
- //
- // An unnecessary empty constructor
- // is added to the class to work
- // around a bug in Compaq's cxx
- // compiler which otherwise reports
- // an error about an omitted
- // initializer for an object of
- // this class further down.
-template <int dim>
-class Beta
-{
- public:
- Beta () {};
- void value_list (const std::vector<Point<dim> > &points,
- std::vector<Point<dim> > &values) const;
-};
+ // @sect4{%Point evaluation}
+
+ // The next thing is to implement
+ // actual evaluation classes. As
+ // noted in the introduction, we'd
+ // like to extract a point value
+ // from the solution, so the first
+ // class does this in its
+ // <code>operator()</code>. The actual point
+ // is given to this class through
+ // the constructor, as well as a
+ // table object into which it will
+ // put its findings.
+ //
+ // Finding out the value of a
+ // finite element field at an
+ // arbitrary point is rather
+ // difficult, if we cannot rely on
+ // knowing the actual finite
+ // element used, since then we
+ // cannot, for example, interpolate
+ // between nodes. For simplicity,
+ // we therefore assume here that
+ // the point at which we want to
+ // evaluate the field is actually a
+ // node. If, in the process of
+ // evaluating the solution, we find
+ // that we did not encounter this
+ // point upon looping over all
+ // vertices, we then have to throw
+ // an exception in order to signal
+ // to the calling functions that
+ // something has gone wrong, rather
+ // than silently ignore this error.
+ //
+ // In the step-9 example program,
+ // we have already seen how such an
+ // exception class can be declared,
+ // using the <code>DeclExceptionN</code>
+ // macros. We use this mechanism
+ // here again.
+ //
+ // From this, the actual
+ // declaration of this class should
+ // be evident. Note that of course
+ // even if we do not list a
+ // destructor explicitely, an
+ // implicit destructor is generated
+ // from the compiler, and it is
+ // virtual just as the one of the
+ // base class.
+ template <int dim>
+ class PointValueEvaluation : public EvaluationBase<dim>
+ {
+ public:
+ PointValueEvaluation (const Point<dim> &evaluation_point,
+ TableHandler &results_table);
+
+ virtual void operator () (const DoFHandler<dim> &dof_handler,
+ const Vector<double> &solution) const;
+
+ DeclException1 (ExcEvaluationPointNotFound,
+ Point<dim>,
+ << "The evaluation point " << arg1
+ << " was not found among the vertices of the present grid.");
+ private:
+ const Point<dim> evaluation_point;
+ TableHandler &results_table;
+ };
+
+
+ // As for the definition, the
+ // constructor is trivial, just
+ // taking data and storing it in
+ // object-local ones:
+ template <int dim>
+ PointValueEvaluation<dim>::
+ PointValueEvaluation (const Point<dim> &evaluation_point,
+ TableHandler &results_table)
+ :
+ evaluation_point (evaluation_point),
+ results_table (results_table)
+ {}
+
- // The implementation of the
- // <code>value_list</code> functions of these
- // classes are rather simple. For
- // simplicity the right hand side is
- // set to be zero but will be
- // assembled anyway.
-template <int dim>
-void RHS<dim>::value_list(const std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int) const
-{
- // Usually we check whether input
- // parameter have the right sizes.
- Assert(values.size()==points.size(),
- ExcDimensionMismatch(values.size(),points.size()));
- for (unsigned int i=0; i<values.size(); ++i)
- values[i]=0;
-}
+ // Now for the function that is
+ // mainly of interest in this
+ // class, the computation of the
+ // point value:
+ template <int dim>
+ void
+ PointValueEvaluation<dim>::
+ operator () (const DoFHandler<dim> &dof_handler,
+ const Vector<double> &solution) const
+ {
+ // First allocate a variable that
+ // will hold the point
+ // value. Initialize it with a
+ // value that is clearly bogus,
+ // so that if we fail to set it
+ // to a reasonable value, we will
+ // note at once. This may not be
+ // necessary in a function as
+ // small as this one, since we
+ // can easily see all possible
+ // paths of execution here, but
+ // it proved to be helpful for
+ // more complex cases, and so we
+ // employ this strategy here as
+ // well.
+ double point_value = 1e20;
+
+ // Then loop over all cells and
+ // all their vertices, and check
+ // whether a vertex matches the
+ // evaluation point. If this is
+ // the case, then extract the
+ // point value, set a flag that
+ // we have found the point of
+ // interest, and exit the loop.
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ bool evaluation_point_found = false;
+ for (; (cell!=endc) && !evaluation_point_found; ++cell)
+ for (unsigned int vertex=0;
+ vertex<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex)
+ if (cell->vertex(vertex) == evaluation_point)
+ {
+ // In order to extract
+ // the point value from
+ // the global solution
+ // vector, pick that
+ // component that belongs
+ // to the vertex of
+ // interest, and, in case
+ // the solution is
+ // vector-valued, take
+ // the first component of
+ // it:
+ point_value = solution(cell->vertex_dof_index(vertex,0));
+ // Note that by this we
+ // have made an
+ // assumption that is not
+ // valid always and
+ // should be documented
+ // in the class
+ // declaration if this
+ // were code for a real
+ // application rather
+ // than a tutorial
+ // program: we assume
+ // that the finite
+ // element used for the
+ // solution we try to
+ // evaluate actually has
+ // degrees of freedom
+ // associated with
+ // vertices. This, for
+ // example, does not hold
+ // for discontinuous
+ // elements, were the
+ // support points for the
+ // shape functions
+ // happen to be located
+ // at the vertices, but
+ // are not associated
+ // with the vertices but
+ // rather with the cell
+ // interior, since
+ // association with
+ // vertices would imply
+ // continuity there. It
+ // would also not hold
+ // for edge oriented
+ // elements, and the
+ // like.
+ //
+ // Ideally, we would
+ // check this at the
+ // beginning of the
+ // function, for example
+ // by a statement like
+ // <code>Assert
+ // (dof_handler.get_fe().dofs_per_vertex
+ // @> 0,
+ // ExcNotImplemented())</code>,
+ // which should make it
+ // quite clear what is
+ // going wrong when the
+ // exception is
+ // triggered. In this
+ // case, we omit it
+ // (which is indeed bad
+ // style), but knowing
+ // that that does not
+ // hurt here, since the
+ // statement
+ // <code>cell-@>vertex_dof_index(vertex,0)</code>
+ // would fail if we asked
+ // it to give us the DoF
+ // index of a vertex if
+ // there were none.
+ //
+ // We stress again that
+ // this restriction on
+ // the allowed finite
+ // elements should be
+ // stated in the class
+ // documentation.
+
+ // Since we found the
+ // right point, we now
+ // set the respective
+ // flag and exit the
+ // innermost loop. The
+ // outer loop will the
+ // also be terminated due
+ // to the set flag.
+ evaluation_point_found = true;
+ break;
+ };
+
+ // Finally, we'd like to make
+ // sure that we have indeed found
+ // the evaluation point, since if
+ // that were not so we could not
+ // give a reasonable value of the
+ // solution there and the rest of
+ // the computations were useless
+ // anyway. So make sure through
+ // the <code>AssertThrow</code> macro
+ // already used in the step-9
+ // program that we have indeed
+ // found this point. If this is
+ // not so, the macro throws an
+ // exception of the type that is
+ // given to it as second
+ // argument, but compared to a
+ // straightforward <code>throw</code>
+ // statement, it fills the
+ // exception object with a set of
+ // additional information, for
+ // example the source file and
+ // line number where the
+ // exception was generated, and
+ // the condition that failed. If
+ // you have a <code>catch</code> clause in
+ // your main function (as this
+ // program has), you will catch
+ // all exceptions that are not
+ // caught somewhere in between
+ // and thus already handled, and
+ // this additional information
+ // will help you find out what
+ // happened and where it went
+ // wrong.
+ AssertThrow (evaluation_point_found,
+ ExcEvaluationPointNotFound(evaluation_point));
+ // Note that we have used the
+ // <code>Assert</code> macro in other
+ // example programs as well. It
+ // differed from the
+ // <code>AssertThrow</code> macro used
+ // here in that it simply aborts
+ // the program, rather than
+ // throwing an exception, and
+ // that it did so only in debug
+ // mode. It was the right macro
+ // to use to check about the size
+ // of vectors passed as arguments
+ // to functions, and the like.
+ //
+ // However, here the situation is
+ // different: whether we find the
+ // evaluation point or not may
+ // change from refinement to
+ // refinement (for example, if
+ // the four cells around point
+ // are coarsened away, then the
+ // point may vanish after
+ // refinement and
+ // coarsening). This is something
+ // that cannot be predicted from
+ // a few number of runs of the
+ // program in debug mode, but
+ // should be checked always, also
+ // in production runs. Thus the
+ // use of the <code>AssertThrow</code>
+ // macro here.
+
+ // Now, if we are sure that we
+ // have found the evaluation
+ // point, we can add the results
+ // into the table of results:
+ results_table.add_value ("DoFs", dof_handler.n_dofs());
+ results_table.add_value ("u(x_0)", point_value);
+ }
- // The flow field is chosen to be
- // circular, counterclockwise, and with
- // the origin as midpoint.
-template <int dim>
-void Beta<dim>::value_list(const std::vector<Point<dim> > &points,
- std::vector<Point<dim> > &values) const
-{
- Assert(values.size()==points.size(),
- ExcDimensionMismatch(values.size(),points.size()));
- for (unsigned int i=0; i<points.size(); ++i)
- {
- const Point<dim> &p=points[i];
- Point<dim> &beta=values[i];
- beta(0) = -p(1);
- beta(1) = p(0);
- beta /= std::sqrt(beta.square());
- }
-}
+ // @sect4{Generating output}
+
+ // A different, maybe slightly odd
+ // kind of <code>evaluation</code> of a
+ // solution is to output it to a
+ // file in a graphical
+ // format. Since in the evaluation
+ // functions we are given a
+ // <code>DoFHandler</code> object and the
+ // solution vector, we have all we
+ // need to do this, so we can do it
+ // in an evaluation class. The
+ // reason for actually doing so
+ // instead of putting it into the
+ // class that computed the solution
+ // is that this way we have more
+ // flexibility: if we choose to
+ // only output certain aspects of
+ // it, or not output it at all. In
+ // any case, we do not need to
+ // modify the solver class, we just
+ // have to modify one of the
+ // modules out of which we build
+ // this program. This form of
+ // encapsulation, as above, helps
+ // us to keep each part of the
+ // program rather simple as the
+ // interfaces are kept simple, and
+ // no access to hidden data is
+ // possible.
+ //
+ // Since this class which generates
+ // the output is derived from the
+ // common <code>EvaluationBase</code> base
+ // class, its main interface is the
+ // <code>operator()</code>
+ // function. Furthermore, it has a
+ // constructor taking a string that
+ // will be used as the base part of
+ // the file name to which output
+ // will be sent (we will augment it
+ // by a number indicating the
+ // number of the refinement cycle
+ // -- the base class has this
+ // information at hand --, and a
+ // suffix), and the constructor
+ // also takes a value that
+ // indicates which format is
+ // requested, i.e. for which
+ // graphics program we shall
+ // generate output (from this we
+ // will then also generate the
+ // suffix of the filename to which
+ // we write).
+ //
+ // Regarding the output format, the
+ // <code>DataOutInterface</code> class
+ // (which is a base class of
+ // <code>DataOut</code> through which we
+ // will access its fields) provides
+ // an enumeration field
+ // <code>OutputFormat</code>, which lists
+ // names for all supported output
+ // formats. At the time of writing
+ // of this program, the supported
+ // graphics formats are represented
+ // by the enum values <code>ucd</code>,
+ // <code>gnuplot</code>, <code>povray</code>,
+ // <code>eps</code>, <code>gmv</code>, <code>tecplot</code>,
+ // <code>tecplot_binary</code>, <code>dx</code>, and
+ // <code>vtk</code>, but this list will
+ // certainly grow over time. Now,
+ // within various functions of that
+ // base class, you can use values
+ // of this type to get information
+ // about these graphics formats
+ // (for example the default suffix
+ // used for files of each format),
+ // and you can call a generic
+ // <code>write</code> function, which then
+ // branches to the
+ // <code>write_gnuplot</code>,
+ // <code>write_ucd</code>, etc functions
+ // which we have used in previous
+ // examples already, based on the
+ // value of a second argument given
+ // to it denoting the required
+ // output format. This mechanism
+ // makes it simple to write an
+ // extensible program that can
+ // decide which output format to
+ // use at runtime, and it also
+ // makes it rather simple to write
+ // the program in a way such that
+ // it takes advantage of newly
+ // implemented output formats,
+ // without the need to change the
+ // application program.
+ //
+ // Of these two fields, the base
+ // name and the output format
+ // descriptor, the constructor
+ // takes values and stores them for
+ // later use by the actual
+ // evaluation function.
+ template <int dim>
+ class SolutionOutput : public EvaluationBase<dim>
+ {
+ public:
+ SolutionOutput (const std::string &output_name_base,
+ const typename DataOut<dim>::OutputFormat output_format);
+
+ virtual void operator () (const DoFHandler<dim> &dof_handler,
+ const Vector<double> &solution) const;
+ private:
+ const std::string output_name_base;
+ const typename DataOut<dim>::OutputFormat output_format;
+ };
+
+
+ template <int dim>
+ SolutionOutput<dim>::
+ SolutionOutput (const std::string &output_name_base,
+ const typename DataOut<dim>::OutputFormat output_format)
+ :
+ output_name_base (output_name_base),
+ output_format (output_format)
+ {}
+
+
+ // After the description above, the
+ // function generating the actual
+ // output is now relatively
+ // straightforward. The only
+ // particularly interesting feature
+ // over previous example programs
+ // is the use of the
+ // <code>DataOut::default_suffix</code>
+ // function, returning the usual
+ // suffix for files of a given
+ // format (e.g. ".eps" for
+ // encapsulated postscript files,
+ // ".gnuplot" for Gnuplot files),
+ // and of the generic
+ // <code>DataOut::write</code> function with
+ // a second argument, which
+ // branches to the actual output
+ // functions for the different
+ // graphics formats, based on the
+ // value of the format descriptor
+ // passed as second argument.
+ //
+ // Also note that we have to prefix
+ // <code>this-@></code> to access a member
+ // variable of the template
+ // dependent base class. The reason
+ // here, and further down in the
+ // program is the same as the one
+ // described in the step-7 example
+ // program (look for <code>two-stage
+ // name lookup</code> there).
+ template <int dim>
+ void
+ SolutionOutput<dim>::operator () (const DoFHandler<dim> &dof_handler,
+ const Vector<double> &solution) const
+ {
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, "solution");
+ data_out.build_patches ();
+
+ std::ostringstream filename;
+ filename << output_name_base << "-"
+ << this->refinement_cycle
+ << data_out.default_suffix (output_format)
+ << std::ends;
+ std::ofstream out (filename.str().c_str());
+
+ data_out.write (out, output_format);
+ }
- // Hence 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) 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.;
- }
+ // @sect4{Other evaluations}
+
+ // In practical applications, one
+ // would add here a list of other
+ // possible evaluation classes,
+ // representing quantities that one
+ // may be interested in. For this
+ // example, that much shall be
+ // sufficient, so we close the
+ // namespace.
}
-
- // @sect3{Class: DGTransportEquation}
+
+ // @sect3{The Laplace solver classes}
+
+ // After defining what we want to
+ // know of the solution, we should
+ // now care how to get at it. We will
+ // pack everything we need into a
+ // namespace of its own, for much the
+ // same reasons as for the
+ // evaluations above.
//
- // Next we define the
- // equation-dependent and
- // DG-method-dependent class
- // <code>DGTransportEquation</code>. Its
- // member functions were already
- // mentioned in the Introduction and
- // will be explained
- // below. Furthermore it includes
- // objects of the previously defined
- // <code>Beta</code>, <code>RHS</code> and
- // <code>BoundaryValues</code> function
- // classes.
-template <int dim>
-class DGTransportEquation
+ // Since we have discussed Laplace
+ // solvers already in considerable
+ // detail in previous examples, there
+ // is not much new stuff
+ // following. Rather, we have to a
+ // great extent cannibalized previous
+ // examples and put them, in slightly
+ // different form, into this example
+ // program. We will therefore mostly
+ // be concerned with discussing the
+ // differences to previous examples.
+ //
+ // Basically, as already said in the
+ // introduction, the lack of new
+ // stuff in this example is
+ // deliberate, as it is more to
+ // demonstrate software design
+ // practices, rather than
+ // mathematics. The emphasis in
+ // explanations below will therefore
+ // be more on the actual
+ // implementation.
+namespace LaplaceSolver
{
- public:
- DGTransportEquation();
-
- void assemble_cell_term(const FEValues<dim>& fe_v,
- FullMatrix<double> &u_v_matrix,
- Vector<double> &cell_vector) const;
-
- void assemble_boundary_term(const FEFaceValues<dim>& fe_v,
- FullMatrix<double> &u_v_matrix,
- Vector<double> &cell_vector) const;
-
- void assemble_face_term1(const FEFaceValuesBase<dim>& fe_v,
- const FEFaceValuesBase<dim>& fe_v_neighbor,
- FullMatrix<double> &u_v_matrix,
- FullMatrix<double> &un_v_matrix) const;
-
- void assemble_face_term2(const FEFaceValuesBase<dim>& fe_v,
- const FEFaceValuesBase<dim>& fe_v_neighbor,
- FullMatrix<double> &u_v_matrix,
- FullMatrix<double> &un_v_matrix,
- FullMatrix<double> &u_vn_matrix,
- FullMatrix<double> &un_vn_matrix) const;
- private:
- const Beta<dim> beta_function;
- const RHS<dim> rhs_function;
- const BoundaryValues<dim> boundary_function;
-};
+ // @sect4{An abstract base class}
+
+ // In defining a Laplace solver, we
+ // start out by declaring an
+ // abstract base class, that has no
+ // functionality itself except for
+ // taking and storing a pointer to
+ // the triangulation to be used
+ // later.
+ //
+ // This base class is very general,
+ // and could as well be used for
+ // any other stationary problem. It
+ // provides declarations of
+ // functions that shall, in derived
+ // classes, solve a problem,
+ // postprocess the solution with a
+ // list of evaluation objects, and
+ // refine the grid,
+ // respectively. None of these
+ // functions actually does
+ // something itself in the base
+ // class.
+ //
+ // Due to the lack of actual
+ // functionality, the programming
+ // style of declaring very abstract
+ // base classes reminds of the
+ // style used in Smalltalk or Java
+ // programs, where all classes are
+ // derived from entirely abstract
+ // classes <code>Object</code>, even number
+ // representations. The author
+ // admits that he does not
+ // particularly like the use of
+ // such a style in C++, as it puts
+ // style over reason. Furthermore,
+ // it promotes the use of virtual
+ // functions for everything (for
+ // example, in Java, all functions
+ // are virtual per se), which,
+ // however, has proven to be rather
+ // inefficient in many applications
+ // where functions are often only
+ // accessing data, not doing
+ // computations, and therefore
+ // quickly return; the overhead of
+ // virtual functions can then be
+ // significant. The opinion of the
+ // author is to have abstract base
+ // classes wherever at least some
+ // part of the code of actual
+ // implementations can be shared
+ // and thus separated into the base
+ // class.
+ //
+ // Besides all these theoretical
+ // questions, we here have a good
+ // reason, which will become
+ // clearer to the reader
+ // below. Basically, we want to be
+ // able to have a family of
+ // different Laplace solvers that
+ // differ so much that no larger
+ // common subset of functionality
+ // could be found. We therefore
+ // just declare such an abstract
+ // base class, taking a pointer to
+ // a triangulation in the
+ // constructor and storing it
+ // henceforth. Since this
+ // triangulation will be used
+ // throughout all computations, we
+ // have to make sure that the
+ // triangulation exists until the
+ // destructor exits. We do this by
+ // keeping a <code>SmartPointer</code> to
+ // this triangulation, which uses a
+ // counter in the triangulation
+ // class to denote the fact that
+ // there is still an object out
+ // there using this triangulation,
+ // thus leading to an abort in case
+ // the triangulation is attempted
+ // to be destructed while this
+ // object still uses it.
+ //
+ // Note that while the pointer
+ // itself is declared constant
+ // (i.e. throughout the lifetime of
+ // this object, the pointer points
+ // to the same object), it is not
+ // declared as a pointer to a
+ // constant triangulation. In fact,
+ // by this we allow that derived
+ // classes refine or coarsen the
+ // triangulation within the
+ // <code>refine_grid</code> function.
+ //
+ // Finally, we have a function
+ // <code>n_dofs</code> is only a tool for
+ // the driver functions to decide
+ // whether we want to go on with
+ // mesh refinement or not. It
+ // returns the number of degrees of
+ // freedom the present simulation
+ // has.
+ template <int dim>
+ class Base
+ {
+ public:
+ Base (Triangulation<dim> &coarse_grid);
+ virtual ~Base ();
+
+ virtual void solve_problem () = 0;
+ virtual void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const = 0;
+ virtual void refine_grid () = 0;
+ virtual unsigned int n_dofs () const = 0;
+
+ protected:
+ const SmartPointer<Triangulation<dim> > triangulation;
+ };
-template <int dim>
-DGTransportEquation<dim>::DGTransportEquation ()
- :
- beta_function (),
- rhs_function (),
- boundary_function ()
-{}
+ // The implementation of the only
+ // two non-abstract functions is
+ // then rather boring:
+ template <int dim>
+ Base<dim>::Base (Triangulation<dim> &coarse_grid)
+ :
+ triangulation (&coarse_grid)
+ {}
- // @sect4{Function: assemble_cell_term}
- //
- // The <code>assemble_cell_term</code>
- // function assembles the cell terms
- // of the discretization.
- // <code>u_v_matrix</code> is a cell matrix,
- // i.e. for a DG method of degree 1,
- // it is of size 4 times 4, and
- // <code>cell_vector</code> is of size 4.
- // When this function is invoked,
- // <code>fe_v</code> is already reinit'ed with the
- // current cell before and includes
- // all shape values needed.
-template <int dim>
-void DGTransportEquation<dim>::assemble_cell_term(
- const FEValues<dim> &fe_v,
- FullMatrix<double> &u_v_matrix,
- Vector<double> &cell_vector) const
-{
- // First we ask <code>fe_v</code> for the
- // shape gradients, shape values and
- // quadrature weights,
- const std::vector<double> &JxW = fe_v.get_JxW_values ();
-
- // Then the flow field beta and the
- // <code>rhs_function</code> are evaluated at
- // the quadrature points,
- std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
- std::vector<double> rhs (fe_v.n_quadrature_points);
-
- beta_function.value_list (fe_v.get_quadrature_points(), beta);
- rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
+ template <int dim>
+ Base<dim>::~Base ()
+ {}
- // and the cell matrix and cell
- // vector are assembled due to the
- // terms $-(u,\beta\cdot\nabla
- // v)_K$ and $(f,v)_K$.
- for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+
+ // @sect4{A general solver class}
+
+ // Following now the main class
+ // that implements assembling the
+ // matrix of the linear system,
+ // solving it, and calling the
+ // postprocessor objects on the
+ // solution. It implements the
+ // <code>solve_problem</code> and
+ // <code>postprocess</code> functions
+ // declared in the base class. It
+ // does not, however, implement the
+ // <code>refine_grid</code> method, as mesh
+ // refinement will be implemented
+ // in a number of derived classes.
+ //
+ // It also declares a new abstract
+ // virtual function,
+ // <code>assemble_rhs</code>, that needs to
+ // be overloaded in subclasses. The
+ // reason is that we will implement
+ // two different classes that will
+ // implement different methods to
+ // assemble the right hand side
+ // vector. This function might also
+ // be interesting in cases where
+ // the right hand side depends not
+ // simply on a continuous function,
+ // but on something else as well,
+ // for example the solution of
+ // another discretized problem,
+ // etc. The latter happens
+ // frequently in non-linear
+ // problems.
+ //
+ // As we mentioned previously, the
+ // actual content of this class is
+ // not new, but a mixture of
+ // various techniques already used
+ // in previous examples. We will
+ // therefore not discuss them in
+ // detail, but refer the reader to
+ // these programs.
+ //
+ // Basically, in a few words, the
+ // constructor of this class takes
+ // pointers to a triangulation, a
+ // finite element, and a function
+ // object representing the boundary
+ // values. These are either passed
+ // down to the base class's
+ // constructor, or are stored and
+ // used to generate a
+ // <code>DoFHandler</code> object
+ // later. Since finite elements and
+ // quadrature formula should match,
+ // it is also passed a quadrature
+ // object.
+ //
+ // The <code>solve_problem</code> sets up
+ // the data structures for the
+ // actual solution, calls the
+ // functions to assemble the linear
+ // system, and solves it.
+ //
+ // The <code>postprocess</code> function
+ // finally takes an evaluation
+ // object and applies it to the
+ // computed solution.
+ //
+ // The <code>n_dofs</code> function finally
+ // implements the pure virtual
+ // function of the base class.
+ template <int dim>
+ class Solver : public virtual Base<dim>
+ {
+ public:
+ Solver (Triangulation<dim> &triangulation,
+ const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const Function<dim> &boundary_values);
+ virtual
+ ~Solver ();
+
+ virtual
+ void
+ solve_problem ();
+
+ virtual
+ void
+ postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
+
+ virtual
+ unsigned int
+ n_dofs () const;
+
+ // In the protected section of
+ // this class, we first have a
+ // number of member variables,
+ // of which the use should be
+ // clear from the previous
+ // examples:
+ protected:
+ const SmartPointer<const FiniteElement<dim> > fe;
+ const SmartPointer<const Quadrature<dim> > quadrature;
+ DoFHandler<dim> dof_handler;
+ Vector<double> solution;
+ const SmartPointer<const Function<dim> > boundary_values;
+
+ // Then we declare an abstract
+ // function that will be used
+ // to assemble the right hand
+ // side. As explained above,
+ // there are various cases for
+ // which this action differs
+ // strongly in what is
+ // necessary, so we defer this
+ // to derived classes:
+ virtual void assemble_rhs (Vector<double> &rhs) const = 0;
+
+ // Next, in the private
+ // section, we have a small
+ // class which represents an
+ // entire linear system, i.e. a
+ // matrix, a right hand side,
+ // and a solution vector, as
+ // well as the constraints that
+ // are applied to it, such as
+ // those due to hanging
+ // nodes. Its constructor
+ // initializes the various
+ // subobjects, and there is a
+ // function that implements a
+ // conjugate gradient method as
+ // solver.
+ private:
+ struct LinearSystem
{
- for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
- u_v_matrix(i,j) -= beta[point]*fe_v.shape_grad(i,point)*
- fe_v.shape_value(j,point) *
- JxW[point];
+ LinearSystem (const DoFHandler<dim> &dof_handler);
+
+ void solve (Vector<double> &solution) const;
- cell_vector(i) += rhs[point] * fe_v.shape_value(i,point) * JxW[point];
- }
-}
+ ConstraintMatrix hanging_node_constraints;
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> matrix;
+ Vector<double> rhs;
+ };
+
+ // Finally, there is a pair of
+ // functions which will be used
+ // to assemble the actual
+ // system matrix. It calls the
+ // virtual function assembling
+ // the right hand side, and
+ // installs a number threads
+ // each running the second
+ // function which assembles
+ // part of the system
+ // matrix. The mechanism for
+ // doing so is the same as in
+ // the step-9 example program.
+ void
+ assemble_linear_system (LinearSystem &linear_system);
+
+ void
+ assemble_matrix (LinearSystem &linear_system,
+ const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
+ const typename DoFHandler<dim>::active_cell_iterator &end_cell,
+ Threads::ThreadMutex &mutex) const;
+ };
+
+
+
+ // Now here comes the constructor
+ // of the class. It does not do
+ // much except store pointers to
+ // the objects given, and generate
+ // <code>DoFHandler</code> object
+ // initialized with the given
+ // pointer to a triangulation. This
+ // causes the DoF handler to store
+ // that pointer, but does not
+ // already generate a finite
+ // element numbering (we only ask
+ // for that in the
+ // <code>solve_problem</code> function).
+ template <int dim>
+ Solver<dim>::Solver (Triangulation<dim> &triangulation,
+ const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const Function<dim> &boundary_values)
+ :
+ Base<dim> (triangulation),
+ fe (&fe),
+ quadrature (&quadrature),
+ dof_handler (triangulation),
+ boundary_values (&boundary_values)
+ {}
+
+
+ // The destructor is simple, it
+ // only clears the information
+ // stored in the DoF handler object
+ // to release the memory.
+ template <int dim>
+ Solver<dim>::~Solver ()
+ {
+ dof_handler.clear ();
+ }
- // @sect4{Function: assemble_boundary_term}
- //
- // The <code>assemble_boundary_term</code>
- // function assembles the face terms
- // at boundary faces. When this
- // function is invoked, <code>fe_v</code> is
- // already reinit'ed with the current
- // cell and current face. Hence it
- // provides the shape values on that
- // boundary face.
-template <int dim>
-void DGTransportEquation<dim>::assemble_boundary_term(
- const FEFaceValues<dim>& fe_v,
- FullMatrix<double> &u_v_matrix,
- Vector<double> &cell_vector) const
-{
- // Again, as in the previous
- // function, we ask the <code>FEValues</code>
- // object for the shape values and
- // the quadrature weights
- const std::vector<double> &JxW = fe_v.get_JxW_values ();
- // but here also for the normals.
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
-
- // We evaluate the flow field
- // and the boundary values at the
- // quadrature points.
- std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
- std::vector<double> g(fe_v.n_quadrature_points);
-
- beta_function.value_list (fe_v.get_quadrature_points(), beta);
- boundary_function.value_list (fe_v.get_quadrature_points(), g);
-
- // Then we assemble cell vector and
- // cell matrix according to the DG
- // method given in the
- // introduction.
- for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
- {
- const double beta_n=beta[point] * normals[point];
- // We assemble the term
- // $(\beta\cdot n
- // u,v)_{\partial K_+}$,
- 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)
- u_v_matrix(i,j) += beta_n *
- fe_v.shape_value(j,point) *
- fe_v.shape_value(i,point) *
- JxW[point];
- else
- // and the term $(\beta\cdot
- // n g,v)_{\partial
- // K_-\cap\partial\Omega}$,
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
- cell_vector(i) -= beta_n *
- g[point] *
- fe_v.shape_value(i,point) *
- JxW[point];
- }
-}
+ // The next function is the one
+ // which delegates the main work in
+ // solving the problem: it sets up
+ // the DoF handler object with the
+ // finite element given to the
+ // constructor of this object, the
+ // creates an object that denotes
+ // the linear system (i.e. the
+ // matrix, the right hand side
+ // vector, and the solution
+ // vector), calls the function to
+ // assemble it, and finally solves
+ // it:
+ template <int dim>
+ void
+ Solver<dim>::solve_problem ()
+ {
+ dof_handler.distribute_dofs (*fe);
+ std::cout << "Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
+
+ solution.reinit (dof_handler.n_dofs());
- // @sect4{Function: assemble_face_term1}
- //
- // The <code>assemble_face_term1</code>
- // function assembles the face terms
- // corresponding to the first version
- // of the DG method, cf. above. For
- // that case, the face terms are
- // given as a sum of integrals over
- // all cell boundaries.
- //
- // When this function is invoked,
- // <code>fe_v</code> and <code>fe_v_neighbor</code> are
- // already reinit'ed with the current
- // cell and the neighoring cell,
- // respectively, as well as with the
- // current face. Hence they provide
- // the inner and outer shape values
- // on the face.
- //
- // In addition to the cell matrix
- // <code>u_v_matrix</code> this function has
- // got a new argument
- // <code>un_v_matrix</code>, that stores
- // contributions to the system matrix
- // that are based on outer values of
- // u, see $\hat u_h$ in the
- // introduction, and inner values of
- // v, see $v_h$. Here we note that
- // <code>un</code> is the short notation for
- // <code>u_neighbor</code> and represents
- // $\hat u_h$.
-template <int dim>
-void DGTransportEquation<dim>::assemble_face_term1(
- const FEFaceValuesBase<dim>& fe_v,
- const FEFaceValuesBase<dim>& fe_v_neighbor,
- FullMatrix<double> &u_v_matrix,
- FullMatrix<double> &un_v_matrix) const
-{
- // Again, as in the previous
- // function, we ask the FEValues
- // objects for the shape values,
- // the quadrature weights and the
- // normals
- const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
-
- // and we evaluate the flow field
- // at the quadrature points.
- std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
- beta_function.value_list (fe_v.get_quadrature_points(), beta);
-
- // Then we assemble the cell
- // matrices according to the DG
- // method given in the
- // introduction.
- for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
- {
- const double beta_n=beta[point] * normals[point];
- // We assemble the term
- // $(\beta\cdot n
- // u,v)_{\partial K_+}$,
- 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)
- u_v_matrix(i,j) += beta_n *
- fe_v.shape_value(j,point) *
- fe_v.shape_value(i,point) *
- JxW[point];
- else
- // and the
- // term $(\beta\cdot n
- // \hat u,v)_{\partial
- // K_-}$.
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
- for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
- un_v_matrix(i,k) += beta_n *
- fe_v_neighbor.shape_value(k,point) *
- fe_v.shape_value(i,point) *
- JxW[point];
- }
-}
+ LinearSystem linear_system (dof_handler);
+ std::cout << "Number of constraints : "
+ << linear_system.hanging_node_constraints.n_constraints()
+ << std::endl;
- // @sect4{Function: assemble_face_term2}
- //
- // Now we look at the
- // <code>assemble_face_term2</code> function
- // that assembles the face terms
- // corresponding to the second
- // version of the DG method,
- // cf. above. For that case the face
- // terms are given as a sum of
- // integrals over all faces. Here we
- // need two additional cell matrices
- // <code>u_vn_matrix</code> and
- // <code>un_vn_matrix</code> that will store
- // contributions due to terms
- // involving u and vn as well as un
- // and vn.
-template <int dim>
-void DGTransportEquation<dim>::assemble_face_term2(
- const FEFaceValuesBase<dim>& fe_v,
- const FEFaceValuesBase<dim>& fe_v_neighbor,
- FullMatrix<double> &u_v_matrix,
- FullMatrix<double> &un_v_matrix,
- FullMatrix<double> &u_vn_matrix,
- FullMatrix<double> &un_vn_matrix) const
-{
- // the first few lines are the same
- const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ assemble_linear_system (linear_system);
+ linear_system.solve (solution);
+ }
- std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
-
- beta_function.value_list (fe_v.get_quadrature_points(), beta);
- for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
- {
- const double beta_n=beta[point] * normals[point];
- if (beta_n>0)
- {
- // This terms 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)
- u_v_matrix(i,j) += beta_n *
- fe_v.shape_value(j,point) *
- fe_v.shape_value(i,point) *
- JxW[point];
-
- // We additionally assemble
- // the term $(\beta\cdot n
- // u,\hat v)_{\partial
- // K_+}$,
- 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)
- u_vn_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)
- un_v_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: $(\beta\cdot n \hat
- // u,\hat v)_{\partial
- // K_-}$.
- 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)
- un_vn_matrix(k,l) -= beta_n *
- fe_v_neighbor.shape_value(l,point) *
- fe_v_neighbor.shape_value(k,point) *
- JxW[point];
- }
- }
-}
+ // As stated above, the
+ // <code>postprocess</code> function takes
+ // an evaluation object, and
+ // applies it to the computed
+ // solution. This function may be
+ // called multiply, once for each
+ // evaluation of the solution which
+ // the user required.
+ template <int dim>
+ void
+ Solver<dim>::
+ postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
+ {
+ postprocessor (dof_handler, solution);
+ }
- // @sect3{Class: DGMethod}
- //
- // After these preparations, we
- // proceed with the main part of this
- // program. The main class, here
- // called <code>DGMethod</code> is basically
- // the main class of step 6. One of
- // the differences is that there's no
- // ConstraintMatrix object. This is,
- // because there are no hanging node
- // constraints in DG discretizations.
-template <int dim>
-class DGMethod
-{
- public:
- DGMethod ();
- ~DGMethod ();
+ // The <code>n_dofs</code> function should
+ // be self-explanatory:
+ template <int dim>
+ unsigned int
+ Solver<dim>::n_dofs () const
+ {
+ return dof_handler.n_dofs();
+ }
+
- void run ();
-
- private:
- void setup_system ();
- void assemble_system1 ();
- void assemble_system2 ();
- void solve (Vector<double> &solution);
- void refine_grid ();
- void output_results (const unsigned int cycle) const;
-
- Triangulation<dim> triangulation;
+ // The following function assembles
+ // matrix and right hand side of
+ // the linear system to be solved
+ // in each step. It goes along the
+ // same lines as used in previous
+ // examples, so we explain it only
+ // briefly:
+ template <int dim>
+ void
+ Solver<dim>::assemble_linear_system (LinearSystem &linear_system)
+ {
+ // First define a convenience
+ // abbreviation for these lengthy
+ // iterator names...
+ typedef
+ typename DoFHandler<dim>::active_cell_iterator
+ active_cell_iterator;
+
+ // ... and use it to split up the
+ // set of cells into a number of
+ // pieces of equal size. The
+ // number of blocks is set to the
+ // default number of threads to
+ // be used, which by default is
+ // set to the number of
+ // processors found in your
+ // computer at startup of the
+ // program:
+ const unsigned int n_threads = multithread_info.n_default_threads;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
+ thread_ranges
+ = Threads::split_range<active_cell_iterator> (dof_handler.begin_active (),
+ dof_handler.end (),
+ n_threads);
+
+ // These ranges are then assigned
+ // to a number of threads which
+ // we create next. Each will
+ // assemble the local cell
+ // matrices on the assigned
+ // cells, and fill the matrix
+ // object with it. Since there is
+ // need for synchronization when
+ // filling the same matrix from
+ // different threads, we need a
+ // mutex here:
+ Threads::ThreadMutex mutex;
+ Threads::ThreadGroup<> threads;
+ for (unsigned int thread=0; thread<n_threads; ++thread)
+ threads += Threads::spawn (*this, &Solver<dim>::assemble_matrix)
+ (linear_system,
+ thread_ranges[thread].first,
+ thread_ranges[thread].second,
+ mutex);
+
+ // While the spawned threads
+ // assemble the system matrix, we
+ // can already compute the right
+ // hand side vector in the main
+ // thread, and condense away the
+ // constraints due to hanging
+ // nodes:
+ assemble_rhs (linear_system.rhs);
+ linear_system.hanging_node_constraints.condense (linear_system.rhs);
+
+ // And while we're already at it
+ // to compute things in parallel,
+ // interpolating boundary values
+ // is one more thing that can be
+ // done independently, so we do
+ // it here:
+ std::map<unsigned int,double> boundary_value_map;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ *boundary_values,
+ boundary_value_map);
- // In contrast to the example code
- // of step-12, this time DG elements
- // of different degree will be used.
- // The different FiniteElement
- // objects for the different polynomial
- // degrees will be stored in the
- // fe_collection object.
- hp::FECollection<dim> fe_collection;
-
- // As already mentioned, the
- // standard DoFHandler has to be
- // replaced by a hp::DoFHandler.
- hp::DoFHandler<dim> dof_handler;
-
- SparsityPattern sparsity;
- SparseMatrix<double> system_matrix;
-
- // We define the quadrature
- // formulae for the cell and the
- // face terms of the
- // discretization.
- // Clearly the hp-Method requires
- // a complete set of quadrature
- // rules for each polynomial
- // degree which will be used in the
- // computations.
- hp::QCollection<dim> quadratures;
- hp::QCollection<dim-1> face_quadratures;
- // And there are two solution
- // vectors, that store the
- // solutions to the problems
- // corresponding to the two
- // different assembling routines
- // <code>assemble_system1</code> and
- // <code>assemble_system2</code>;
- Vector<double> solution1;
- Vector<double> solution2;
- Vector<double> right_hand_side;
+ // If this is done, wait for the
+ // matrix assembling threads, and
+ // condense the constraints in
+ // the matrix as well:
+ threads.join_all ();
+ linear_system.hanging_node_constraints.condense (linear_system.matrix);
+
+ // Now that we have the linear
+ // system, we can also treat
+ // boundary values, which need to
+ // be eliminated from both the
+ // matrix and the right hand
+ // side:
+ MatrixTools::apply_boundary_values (boundary_value_map,
+ linear_system.matrix,
+ solution,
+ linear_system.rhs);
+
+ }
+
+
+ // The second of this pair of
+ // functions takes a range of cell
+ // iterators, and assembles the
+ // system matrix on this part of
+ // the domain. Since it's actions
+ // have all been explained in
+ // previous programs, we do not
+ // comment on it any more, except
+ // for one pointe below.
+ template <int dim>
+ void
+ Solver<dim>::assemble_matrix (LinearSystem &linear_system,
+ const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
+ const typename DoFHandler<dim>::active_cell_iterator &end_cell,
+ Threads::ThreadMutex &mutex) const
+ {
+ FEValues<dim> fe_values (*fe, *quadrature,
+ update_gradients | update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe->dofs_per_cell;
+ const unsigned int n_q_points = quadrature->n_quadrature_points;
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ for (typename DoFHandler<dim>::active_cell_iterator cell=begin_cell;
+ cell!=end_cell; ++cell)
+ {
+ cell_matrix = 0;
+
+ fe_values.reinit (cell);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point) *
+ fe_values.JxW(q_point));
+
+
+ cell->get_dof_indices (local_dof_indices);
+
+ // In the step-9 program, we
+ // have shown that you have
+ // to use the mutex to lock
+ // the matrix when copying
+ // the elements from the
+ // local to the global
+ // matrix. This was necessary
+ // to avoid that two threads
+ // access it at the same
+ // time, eventually
+ // overwriting their
+ // respective
+ // work. Previously, we have
+ // used the <code>acquire</code> and
+ // <code>release</code> functions of
+ // the mutex to lock and
+ // unlock the mutex,
+ // respectively. While this
+ // is valid, there is one
+ // possible catch: if between
+ // the locking operation and
+ // the unlocking operation an
+ // exception is thrown, the
+ // mutex remains in the
+ // locked state, and in some
+ // cases this might lead to
+ // deadlocks. A similar
+ // situation arises, when one
+ // changes the code to have a
+ // return statement somewhere
+ // in the middle of the
+ // locked block, and forgets
+ // that before we call
+ // <code>return</code>, we also have
+ // to unlock the mutex. This
+ // all is not be a problem
+ // here, but we want to show
+ // the general technique to
+ // cope with these problems
+ // nevertheless: have an
+ // object that upon
+ // initialization (i.e. in
+ // its constructor) locks the
+ // mutex, and on running the
+ // destructor unlocks it
+ // again. This is called the
+ // <code>scoped lock</code> pattern
+ // (apparently invented by
+ // Doug Schmidt originally),
+ // and it works because
+ // destructors of local
+ // objects are also run when
+ // we exit the function
+ // either through a
+ // <code>return</code> statement, or
+ // when an exception is
+ // raised. Thus, it is
+ // guaranteed that the mutex
+ // will always be unlocked
+ // when we exit this part of
+ // the program, whether the
+ // operation completed
+ // successfully or not,
+ // whether the exit path was
+ // something we implemented
+ // willfully or whether the
+ // function was exited by an
+ // exception that we did not
+ // forsee.
+ //
+ // deal.II implements the
+ // scoped locking pattern in
+ // the
+ // ThreadMutex::ScopedLock
+ // class: it takes the mutex
+ // in the constructor and
+ // locks it; in its
+ // destructor, it unlocks it
+ // again. So here is how it
+ // is used:
+ Threads::ThreadMutex::ScopedLock lock (mutex);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ linear_system.matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+ // Here, at the brace, the
+ // current scope ends, so the
+ // <code>lock</code> variable goes out
+ // of existence and its
+ // destructor the mutex is
+ // unlocked.
+ };
+ }
+
+
+ // Now for the functions that
+ // implement actions in the linear
+ // system class. First, the
+ // constructor initializes all data
+ // elements to their correct sizes,
+ // and sets up a number of
+ // additional data structures, such
+ // as constraints due to hanging
+ // nodes. Since setting up the
+ // hanging nodes and finding out
+ // about the nonzero elements of
+ // the matrix is independent, we do
+ // that in parallel (if the library
+ // was configured to use
+ // concurrency, at least;
+ // otherwise, the actions are
+ // performed sequentially). Note
+ // that we spawn only one thread,
+ // and do the second action in the
+ // main thread. Since only one
+ // thread is generated, we don't
+ // use the <code>Threads::ThreadGroup</code>
+ // class here, but rather use the
+ // one created thread object
+ // directly to wait for this
+ // particular thread's exit.
+ //
+ // Note that taking up the address
+ // of the
+ // <code>DoFTools::make_hanging_node_constraints</code>
+ // function is a little tricky,
+ // since there are actually three
+ // of them, one for each supported
+ // space dimension. Taking
+ // addresses of overloaded
+ // functions is somewhat
+ // complicated in C++, since the
+ // address-of operator <code>&</code> in
+ // that case returns more like a
+ // set of values (the addresses of
+ // all functions with that name),
+ // and selecting the right one is
+ // then the next step. If the
+ // context dictates which one to
+ // take (for example by assigning
+ // to a function pointer of known
+ // type), then the compiler can do
+ // that by itself, but if this set
+ // of pointers shall be given as
+ // the argument to a function that
+ // takes a template, the compiler
+ // could choose all without having
+ // a preference for one. We
+ // therefore have to make it clear
+ // to the compiler which one we
+ // would like to have; for this, we
+ // could use a cast, but for more
+ // clarity, we assign it to a
+ // temporary <code>mhnc_p</code> (short for
+ // <code>pointer to
+ // make_hanging_node_constraints</code>)
+ // with the right type, and using
+ // this pointer instead.
+ template <int dim>
+ Solver<dim>::LinearSystem::
+ LinearSystem (const DoFHandler<dim> &dof_handler)
+ {
+ hanging_node_constraints.clear ();
+
+ void (*mhnc_p) (const DoFHandler<dim> &,
+ ConstraintMatrix &)
+ = &DoFTools::make_hanging_node_constraints;
- // Finally this class includes an
- // object of the
- // DGTransportEquations class
- // described above.
- const DGTransportEquation<dim> dg;
-};
+ Threads::Thread<>
+ mhnc_thread = Threads::spawn (mhnc_p)(dof_handler,
+ hanging_node_constraints);
+
+ sparsity_pattern.reinit (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+ // Wait until the
+ // <code>hanging_node_constraints</code>
+ // object is fully set up, then
+ // close it and use it to
+ // condense the sparsity pattern:
+ mhnc_thread.join ();
+ hanging_node_constraints.close ();
+ hanging_node_constraints.condense (sparsity_pattern);
+
+ // Finally, close the sparsity
+ // pattern, initialize the
+ // matrix, and set the right hand
+ // side vector to the right size.
+ sparsity_pattern.compress();
+ matrix.reinit (sparsity_pattern);
+ rhs.reinit (dof_handler.n_dofs());
+ }
-template <int dim>
-DGMethod<dim>::DGMethod ()
- :
- dof_handler (triangulation),
- dg ()
-{
- // Change here for hp
- // methods of
- // different maximum degrees.
- const unsigned int hp_degree = 5;
- for (unsigned int i = 0; i < hp_degree; ++i)
- {
- fe_collection.push_back (FE_DGQ<dim> (i));
- quadratures.push_back (QGauss<dim> (i+2));
- face_quadratures.push_back (QGauss<dim-1> (i+2));
- }
-}
+ // The second function of this
+ // class simply solves the linear
+ // system by a preconditioned
+ // conjugate gradient method. This
+ // has been extensively discussed
+ // before, so we don't dwell into
+ // it any more.
+ template <int dim>
+ void
+ Solver<dim>::LinearSystem::solve (Vector<double> &solution) const
+ {
+ SolverControl solver_control (1000, 1e-12);
+ SolverCG<> cg (solver_control);
-template <int dim>
-DGMethod<dim>::~DGMethod ()
-{
- dof_handler.clear ();
-}
+ PreconditionSSOR<> preconditioner;
+ preconditioner.initialize(matrix, 1.2);
+ cg.solve (matrix, solution, rhs, preconditioner);
-template <int dim>
-void DGMethod<dim>::setup_system ()
-{
- // First we need to distribute the
- // DoFs.
- dof_handler.distribute_dofs (fe_collection);
- // In order to get a good
- // preconditioner, the degrees of
- // freedom should be ordered in
- // downstream direction. First, we
- // initalize a vector fairly close
- // to the real vector field; since
- // this is for preconditioning
- // only, a rough approximation is
- // sufficient.
- Point<dim> sorting_direction;
- for (unsigned int d=0;d<dim;++d)
- sorting_direction(d) = 1.;
- // Now do the sorting of the
- // degrees of freedom.
- DoFRenumbering::downstream_dg(dof_handler, sorting_direction);
-
- // The DoFs of a cell are coupled
- // with all DoFs of all neighboring
- // cells. Therefore the maximum
- // number of matrix entries per row
- // is needed when all neighbors of
- // a cell are once more refined
- // than the cell under
- // consideration.
- CompressedSparsityPattern compressed_pattern (dof_handler.n_dofs ());
- DoFTools::make_sparsity_pattern (dof_handler, compressed_pattern);
- DoFTools::make_flux_sparsity_pattern (dof_handler, compressed_pattern);
-
- sparsity.copy_from(compressed_pattern);
- system_matrix.reinit (sparsity);
-
- solution1.reinit (dof_handler.n_dofs());
- solution2.reinit (dof_handler.n_dofs());
- right_hand_side.reinit (dof_handler.n_dofs());
-}
+ hanging_node_constraints.distribute (solution);
+ }
- // @sect4{Function: assemble_system1}
- //
- // We proceed with the
- // <code>assemble_system1</code> function that
- // implements the DG discretization
- // in its first version. This
- // function repeatedly calls the
- // <code>assemble_cell_term</code>,
- // <code>assemble_boundary_term</code> and
- // <code>assemble_face_term1</code> functions
- // of the <code>DGTransportEquation</code>
- // object. The
- // <code>assemble_boundary_term</code> covers
- // the first case mentioned in the
- // introduction.
- //
- // 1. face is at boundary
- //
- // This function takes a
- // <code>FEFaceValues</code> object as
- // argument. In contrast to that
- // <code>assemble_face_term1</code>
- // takes two <code>FEFaceValuesBase</code>
- // objects; one for the shape
- // functions on the current cell and
- // the other for shape functions on
- // the neighboring cell under
- // consideration. Both objects are
- // either of class <code>FEFaceValues</code>
- // or of class <code>FESubfaceValues</code>
- // (both derived from
- // <code>FEFaceValuesBase</code>) according to
- // the remaining cases mentioned
- // in the introduction:
- //
- // 2. neighboring cell is finer
- // (current cell: <code>FESubfaceValues</code>,
- // neighboring cell: <code>FEFaceValues</code>);
- //
- // 3. neighboring cell is of the same
- // refinement level (both, current
- // and neighboring cell:
- // <code>FEFaceValues</code>);
- //
- // 4. neighboring cell is coarser
- // (current cell: <code>FEFaceValues</code>,
- // neighboring cell:
- // <code>FESubfaceValues</code>).
- //
- // If we considered globally refined
- // meshes then only case 3 would
- // occur. But as we consider also
- // locally refined meshes we need to
- // distinguish all four cases making
- // the following assembling function
- // a bit longish.
-template <int dim>
-void DGMethod<dim>::assemble_system1 ()
-{
- // First we create the
- // <code>UpdateFlags</code> for the
- // <code>FEValues</code> and the
- // <code>FEFaceValues</code> objects.
- const UpdateFlags update_flags = update_values
- | update_gradients
- | update_q_points
- | update_JxW_values;
-
- // Note, that on faces we do not
- // need gradients but we need
- // normal vectors.
- const UpdateFlags face_update_flags = update_values
- | update_q_points
- | update_JxW_values
- | update_normal_vectors;
-
- // On the neighboring cell we only
- // need the shape values. Given a
- // specific face, the quadrature
- // points and `JxW values' are the
- // same as for the current cells,
- // the normal vectors are known to
- // be the negative of the normal
- // vectors of the current cell.
- const UpdateFlags neighbor_face_update_flags = update_values;
-
- // Then we create the <code>FEValues</code>
- // object. Here, we use the default
- // MappingQ1. different mapping
- // create a MappingCollection first
- // and call the respective
- // hp::FEValues constructor.
- hp::FEValues<dim> fe_v_x (fe_collection, quadratures, update_flags);
-
- // Similarly we create the
- // <code>FEFaceValues</code> and
- // <code>FESubfaceValues</code> objects for
- // both, the current and the
- // neighboring cell. Within the
- // following nested loop over all
- // cells and all faces of the cell
- // they will be reinited to the
- // current cell and the face (and
- // subface) number.
- hp::FEFaceValues<dim> fe_v_face_x (
- fe_collection, face_quadratures, face_update_flags);
- hp::FESubfaceValues<dim> fe_v_subface_x (
- fe_collection, face_quadratures, face_update_flags);
- hp::FEFaceValues<dim> fe_v_face_neighbor_x (
- fe_collection, face_quadratures, neighbor_face_update_flags);
- hp::FESubfaceValues<dim> fe_v_subface_neighbor_x (
- fe_collection, face_quadratures, neighbor_face_update_flags);
-
- // Now we create the cell matrices
- // and vectors. Here we need two
- // cell matrices, both for face
- // terms that include test
- // functions <code>v</code> (shape functions
- // of the current cell). To be more
- // precise, the first matrix will
- // include the `u and v terms' and
- // the second that will include the
- // `un and v terms'. Here we recall
- // the convention that `un' is
- // the shortcut for `u_neighbor'
- // and represents the $u_hat$, see
- // introduction.
- const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell ();
-
- FullMatrix<double> u_v_matrix (max_dofs_per_cell, max_dofs_per_cell);
- FullMatrix<double> un_v_matrix (max_dofs_per_cell, max_dofs_per_cell);
- Vector<double> cell_vector (max_dofs_per_cell);
-
- // Furthermore we need some cell
- // iterators.
- typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
-
- // Now we start the loop over all
- // active cells.
- for (;cell!=endc; ++cell)
- {
- const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
- std::vector<unsigned int> dofs (dofs_per_cell);
- std::vector<unsigned int> dofs_neighbor;
-
- // In the
- // <code>assemble_face_term1</code>
- // function contributions to
- // the cell matrices and the
- // cell vector are only
- // ADDED. Therefore on each
- // cell we need to reset the
- // <code>u_v_matrix</code> and
- // <code>cell_vector</code> to zero,
- // before assembling the cell terms.
- u_v_matrix = 0;
- cell_vector = 0;
-
- // Now we reinit the <code>FEValues</code>
- // object for the current cell
- fe_v_x.reinit (cell);
-
- // and call the function
- // that assembles the cell
- // terms. The first argument is
- // the <code>FEValues</code> that was
- // previously reinit'ed on the
- // current cell.
- dg.assemble_cell_term(fe_v_x.get_present_fe_values (),
- u_v_matrix,
- cell_vector);
-
- // As in previous examples the
- // vector `dofs' includes the
- // dof_indices.
- dofs.resize (dofs_per_cell);
- cell->get_dof_indices (dofs);
-
- // This is the start of the
- // nested loop over all faces.
- for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
- {
- // First we set the face
- // iterator
- typename hp::DoFHandler<dim>::face_iterator face=cell->face(face_no);
-
- // and clear the
- // <code>un_v_matrix</code> on each
- // face.
- un_v_matrix = 0;
-
- // Now we distinguish the
- // four different cases in
- // the ordering mentioned
- // above. We start with
- // faces belonging to the
- // boundary of the domain.
- if (face->at_boundary())
- {
- // We reinit the
- // <code>FEFaceValues</code>
- // object to the
- // current face
- fe_v_face_x.reinit (cell, face_no);
-
- // and assemble the
- // corresponding face
- // terms.
- dg.assemble_boundary_term(fe_v_face_x.get_present_fe_values (),
- u_v_matrix,
- cell_vector);
- }
- else
- {
- // Now we are not on
- // the boundary of the
- // domain, therefore
- // there must exist a
- // neighboring cell.
- typename hp::DoFHandler<dim>::cell_iterator neighbor=
- cell->neighbor(face_no);
-
- // We proceed with the
- // second and most
- // complicated case:
- // the neighboring cell
- // is more refined than
- // the current cell. As
- // in deal.II
- // neighboring cells
- // are restricted to
- // have a level
- // difference of not
- // more than one, the
- // neighboring cell is
- // known to be at most
- // ONCE more refined
- // than the current
- // cell. Furthermore
- // also the face is
- // more refined,
- // i.e. it has
- // children. Here we
- // note that the
- // following part of
- // code will not work
- // for <code>dim==1</code>.
- if (face->has_children())
- {
- // First we store
- // which number the
- // current cell has
- // in the list of
- // neighbors of the
- // neighboring
- // cell. Hence,
- // neighbor-@>neighbor(neighbor2)
- // equals the
- // current cell
- // <code>cell</code>.
- const unsigned int neighbor2=
- cell->neighbor_of_neighbor(face_no);
-
-
- // We loop over
- // subfaces
- for (unsigned int subface_no=0;
- subface_no<face->n_children(); ++subface_no)
- {
- // and set the
- // cell
- // iterator
- // <code>neighbor_child</code>
- // to the cell
- // placed
- // `behind' the
- // current
- // subface.
- typename hp::DoFHandler<dim>::active_cell_iterator neighbor_child=
- neighbor->child(GeometryInfo<dim>::
- child_cell_on_face(neighbor2,subface_no));
-
- // an additional speciality
- // for the hp method appears
- // on the faces. To get an
- // efficient assembly, the
- // lowest order but
- // sufficient quadrature
- // rule should be used. Hence
- // the face quadrature rule of the
- // higher order element
- // will be used.
- const unsigned int quadrature_index =
- std::max (neighbor_child->active_fe_index (),
- cell->active_fe_index ());
-
-
- // As these are
- // quite
- // complicated
- // indirections
- // which one
- // does not
- // usually get
- // right at
- // first
- // attempt we
- // check for
- // the internal
- // consistency.
- Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
- ExcInternalError());
- Assert (!neighbor_child->has_children(), ExcInternalError());
-
- // We need to
- // reset the
- // <code>un_v_matrix</code>
- // on each
- // subface
- // because on
- // each subface
- // the <code>un</code>
- // belong to
- // different
- // neighboring
- // cells.
- un_v_matrix = 0;
-
- // As already
- // mentioned
- // above for
- // the current
- // case (case
- // 2) we employ
- // the
- // <code>FESubfaceValues</code>
- // of the
- // current
- // cell (here
- // reinited for
- // the current
- // cell, face
- // and subface)
- // and we
- // employ the
- // FEFaceValues
- // of the
- // neighboring
- // child cell.
- fe_v_subface_x.reinit (cell, face_no, subface_no, quadrature_index);
- fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, quadrature_index);
-
- dg.assemble_face_term1(fe_v_subface_x.get_present_fe_values (),
- fe_v_face_neighbor_x.get_present_fe_values (),
- u_v_matrix,
- un_v_matrix);
-
- // Then we get
- // the dof
- // indices of
- // the
- // neighbor_child
- // cell
- dofs_neighbor.resize (neighbor_child->get_fe().dofs_per_cell);
- neighbor_child->get_dof_indices (dofs_neighbor);
-
- // and
- // distribute
- // <code>un_v_matrix</code>
- // to the
- // system_matrix
- for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- for (unsigned int k=0; k<neighbor_child->get_fe().dofs_per_cell; ++k)
- system_matrix.add(dofs[i], dofs_neighbor[k],
- un_v_matrix(i,k));
- }
- // End of <code>if
- // (face-@>has_children())</code>
- }
- else
- {
- // We proceed with
- // case 3,
- // i.e. neighboring
- // cell is of the
- // same refinement
- // level as the
- // current cell.
- if (neighbor->level() == cell->level())
- {
- // Like before
- // we store
- // which number
- // the current
- // cell has in
- // the list of
- // neighbors of
- // the
- // neighboring
- // cell.
- const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
-
- // Like before. Use
- // quadrature rule
- // of higher order
- // cell.
- const unsigned int quadrature_index =
- std::max (neighbor->active_fe_index (),
- cell->active_fe_index ());
-
- // We reinit
- // the
- // <code>FEFaceValues</code>
- // of the
- // current and
- // neighboring
- // cell to the
- // current face
- // and assemble
- // the
- // corresponding
- // face terms.
- fe_v_face_x.reinit (cell, face_no, quadrature_index);
- fe_v_face_neighbor_x.reinit (neighbor, neighbor2, quadrature_index);
-
- dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (),
- fe_v_face_neighbor_x.get_present_fe_values (),
- u_v_matrix,
- un_v_matrix);
- // End of <code>if
- // (neighbor-@>level()
- // ==
- // cell-@>level())</code>
- }
- else
- {
- // Finally we
- // consider
- // case 4. When
- // the
- // neighboring
- // cell is not
- // finer and
- // not of the
- // same
- // refinement
- // level as the
- // current cell
- // it must be
- // coarser.
- Assert(neighbor->level() < cell->level(), ExcInternalError());
-
- // Find out the
- // how many'th
- // face_no and
- // subface_no
- // the current
- // face is
- // w.r.t. the
- // neighboring
- // cell.
- const std::pair<unsigned int, unsigned int> faceno_subfaceno=
- cell->neighbor_of_coarser_neighbor(face_no);
- const unsigned int neighbor_face_no=faceno_subfaceno.first,
- neighbor_subface_no=faceno_subfaceno.second;
-
- Assert (neighbor->neighbor(neighbor_face_no)
- ->child(GeometryInfo<dim>::child_cell_on_face(
- face_no,neighbor_subface_no)) == cell, ExcInternalError());
-
-
- // Like before. Use
- // quadrature rule
- // of higher order
- // cell.
- const unsigned int quadrature_index =
- std::max (neighbor->active_fe_index (),
- cell->active_fe_index ());
-
- // Reinit the
- // appropriate
- // <code>FEFaceValues</code>
- // and assemble
- // the face
- // terms.
- fe_v_face_x.reinit (cell, face_no, quadrature_index);
- fe_v_subface_neighbor_x.reinit (neighbor, neighbor_face_no,
- neighbor_subface_no, quadrature_index);
-
- dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (),
- fe_v_subface_neighbor_x.get_present_fe_values (),
- u_v_matrix,
- un_v_matrix);
- }
-
- // Now we get the
- // dof indices of
- // the
- // <code>neighbor_child</code>
- // cell,
- dofs_neighbor.resize (neighbor->get_fe().dofs_per_cell);
- neighbor->get_dof_indices (dofs_neighbor);
-
- // and distribute the
- // <code>un_v_matrix</code>.
- for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- for (unsigned int k=0; k<neighbor->get_fe().dofs_per_cell; ++k)
- system_matrix.add(dofs[i], dofs_neighbor[k],
- un_v_matrix(i,k));
- }
- // End of <code>face not at boundary</code>:
- }
- // End of loop over all faces:
- }
-
- // Finally we distribute the
- // <code>u_v_matrix</code>
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- system_matrix.add(dofs[i], dofs[j], u_v_matrix(i,j));
+
+
+ // @sect4{A primal solver}
+
+ // In the previous section, a base
+ // class for Laplace solvers was
+ // implemented, that lacked the
+ // functionality to assemble the
+ // right hand side vector, however,
+ // for reasons that were explained
+ // there. Now we implement a
+ // corresponding class that can do
+ // this for the case that the right
+ // hand side of a problem is given
+ // as a function object.
+ //
+ // The actions of the class are
+ // rather what you have seen
+ // already in previous examples
+ // already, so a brief explanation
+ // should suffice: the constructor
+ // takes the same data as does that
+ // of the underlying class (to
+ // which it passes all information)
+ // except for one function object
+ // that denotes the right hand side
+ // of the problem. A pointer to
+ // this object is stored (again as
+ // a <code>SmartPointer</code>, in order to
+ // make sure that the function
+ // object is not deleted as long as
+ // it is still used by this class).
+ //
+ // The only functional part of this
+ // class is the <code>assemble_rhs</code>
+ // method that does what its name
+ // suggests.
+ template <int dim>
+ class PrimalSolver : public Solver<dim>
+ {
+ public:
+ PrimalSolver (Triangulation<dim> &triangulation,
+ const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const Function<dim> &rhs_function,
+ const Function<dim> &boundary_values);
+ protected:
+ const SmartPointer<const Function<dim> > rhs_function;
+ virtual void assemble_rhs (Vector<double> &rhs) const;
+ };
+
+
+ // The constructor of this class
+ // basically does what it is
+ // announced to do above...
+ template <int dim>
+ PrimalSolver<dim>::
+ PrimalSolver (Triangulation<dim> &triangulation,
+ const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const Function<dim> &rhs_function,
+ const Function<dim> &boundary_values)
+ :
+ Base<dim> (triangulation),
+ Solver<dim> (triangulation, fe,
+ quadrature, boundary_values),
+ rhs_function (&rhs_function)
+ {}
+
+
+
+ // ... as does the <code>assemble_rhs</code>
+ // function. Since this is
+ // explained in several of the
+ // previous example programs, we
+ // leave it at that.
+ template <int dim>
+ void
+ PrimalSolver<dim>::
+ assemble_rhs (Vector<double> &rhs) const
+ {
+ FEValues<dim> fe_values (*this->fe, *this->quadrature,
+ update_values | update_q_points |
+ update_JxW_values);
+
+ const unsigned int dofs_per_cell = this->fe->dofs_per_cell;
+ const unsigned int n_q_points = this->quadrature->n_quadrature_points;
+
+ Vector<double> cell_rhs (dofs_per_cell);
+ std::vector<double> rhs_values (n_q_points);
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = this->dof_handler.begin_active(),
+ endc = this->dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ cell_rhs = 0;
+ fe_values.reinit (cell);
+ rhs_function->value_list (fe_values.get_quadrature_points(),
+ rhs_values);
- // and the cell vector.
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- right_hand_side(dofs[i]) += cell_vector(i);
- }
-}
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ rhs_values[q_point] *
+ fe_values.JxW(q_point));
+
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ rhs(local_dof_indices[i]) += cell_rhs(i);
+ };
+ }
- // @sect4{Function: assemble_system2}
- //
- // We proceed with the
- // <code>assemble_system2</code> function that
- // implements the DG discretization
- // in its second version. This
- // function is very similar to the
- // <code>assemble_system1</code>
- // function. Therefore, here we only
- // discuss the differences between
- // the two functions. This function
- // repeatedly calls the
- // <code>assemble_face_term2</code> function
- // of the DGTransportEquation object,
- // that assembles the face terms
- // written as a sum of integrals over
- // all faces. Therefore, we need to
- // make sure that each face is
- // treated only once. This is achieved
- // by introducing the rule:
- //
- // a) If the current and the
- // neighboring cells are of the same
- // refinement level we access and
- // treat the face from the cell with
- // lower index.
- //
- // b) If the two cells are of
- // different refinement levels we
- // access and treat the face from the
- // coarser cell.
- //
- // Due to rule b) we do not need to
- // consider case 4 (neighboring cell
- // is coarser) any more.
+ // @sect4{Local refinement by the Kelly error indicator}
+
+ // The second class implementing
+ // refinement strategies uses the
+ // Kelly refinemet indicator used
+ // in various example programs
+ // before. Since this indicator is
+ // already implemented in a class
+ // of its own inside the deal.II
+ // library, there is not much t do
+ // here except cal the function
+ // computing the indicator, then
+ // using it to select a number of
+ // cells for refinement and
+ // coarsening, and refinement the
+ // mesh accordingly.
+ //
+ // Again, this should now be
+ // sufficiently standard to allow
+ // the omission of further
+ // comments.
+ template <int dim>
+ class RefinementKelly : public PrimalSolver<dim>
+ {
+ public:
+ RefinementKelly (Triangulation<dim> &coarse_grid,
+ const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const Function<dim> &rhs_function,
+ const Function<dim> &boundary_values);
-template <int dim>
-void DGMethod<dim>::assemble_system2 ()
-{
- const UpdateFlags update_flags = update_values
- | update_gradients
- | update_q_points
- | update_JxW_values;
-
- const UpdateFlags face_update_flags = update_values
- | update_q_points
- | update_JxW_values
- | update_normal_vectors;
-
- const UpdateFlags neighbor_face_update_flags = update_values;
-
- // Here we do not need
- // <code>fe_v_face_neighbor</code> as case 4
- // does not occur.
- hp::FEValues<dim> fe_v_x (
- fe_collection, quadratures, update_flags);
- hp::FEFaceValues<dim> fe_v_face_x (
- fe_collection, face_quadratures, face_update_flags);
- hp::FESubfaceValues<dim> fe_v_subface_x (
- fe_collection, face_quadratures, face_update_flags);
- hp::FEFaceValues<dim> fe_v_face_neighbor_x (
- fe_collection, face_quadratures, neighbor_face_update_flags);
-
- const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell ();
-
- FullMatrix<double> u_v_matrix (max_dofs_per_cell, max_dofs_per_cell);
- FullMatrix<double> un_v_matrix (max_dofs_per_cell, max_dofs_per_cell);
-
- // Additionally we need the
- // following two cell matrices,
- // both for face term that include
- // test function <code>vn</code> (shape
- // functions of the neighboring
- // cell). To be more precise, the
- // first matrix will include the `u
- // and vn terms' and the second
- // that will include the `un and vn
- // terms'.
- FullMatrix<double> u_vn_matrix (max_dofs_per_cell, max_dofs_per_cell);
- FullMatrix<double> un_vn_matrix (max_dofs_per_cell, max_dofs_per_cell);
-
- Vector<double> cell_vector (max_dofs_per_cell);
-
- // The following lines are roughly
- // the same as in the previous
- // function.
- typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (;cell!=endc; ++cell)
- {
- const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
- std::vector<unsigned int> dofs (dofs_per_cell);
- std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
+ virtual void refine_grid ();
+ };
- u_v_matrix = 0;
- cell_vector = 0;
- fe_v_x.reinit (cell);
- dg.assemble_cell_term(fe_v_x.get_present_fe_values (),
- u_v_matrix,
- cell_vector);
-
- cell->get_dof_indices (dofs);
+ template <int dim>
+ RefinementKelly<dim>::
+ RefinementKelly (Triangulation<dim> &coarse_grid,
+ const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const Function<dim> &rhs_function,
+ const Function<dim> &boundary_values)
+ :
+ Base<dim> (coarse_grid),
+ PrimalSolver<dim> (coarse_grid, fe, quadrature,
+ rhs_function, boundary_values)
+ {}
+
+
+
+ template <int dim>
+ void
+ RefinementKelly<dim>::refine_grid ()
+ {
+ Vector<float> estimated_error_per_cell (this->triangulation->n_active_cells());
+ KellyErrorEstimator<dim>::estimate (this->dof_handler,
+ QGauss<dim-1>(3),
+ typename FunctionMap<dim>::type(),
+ this->solution,
+ estimated_error_per_cell);
+ GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.03);
+ this->triangulation->execute_coarsening_and_refinement ();
+ }
- for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
- {
- typename hp::DoFHandler<dim>::face_iterator face=
- cell->face(face_no);
-
- // Case 1:
- if (face->at_boundary())
- {
- fe_v_face_x.reinit (cell, face_no);
-
- dg.assemble_boundary_term(fe_v_face_x.get_present_fe_values (),
- u_v_matrix,
- cell_vector);
- }
- else
- {
- Assert (cell->neighbor(face_no).state() == IteratorState::valid,
- ExcInternalError());
- typename hp::DoFHandler<dim>::cell_iterator neighbor=
- cell->neighbor(face_no);
-
- const unsigned int dofs_on_neighbor = neighbor->get_fe().dofs_per_cell;
-
- // Case 2:
- if (face->has_children())
- {
- const unsigned int neighbor2=
- cell->neighbor_of_neighbor(face_no);
-
- for (unsigned int subface_no=0;
- subface_no<face->n_children(); ++subface_no)
- {
- typename hp::DoFHandler<dim>::cell_iterator neighbor_child=
- neighbor->child(GeometryInfo<dim>::child_cell_on_face(
- neighbor2,subface_no));
- const unsigned int dofs_on_neighbor_child = neighbor_child->get_fe().dofs_per_cell;
-
- const unsigned int quadrature_index =
- std::max (neighbor_child->active_fe_index (),
- cell->active_fe_index ());
-
- Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
- ExcInternalError());
- Assert (!neighbor_child->has_children(), ExcInternalError());
-
- un_v_matrix = 0;
- u_vn_matrix = 0;
- un_vn_matrix = 0;
-
- fe_v_subface_x.reinit (cell, face_no, subface_no, quadrature_index);
- fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, quadrature_index);
-
- dg.assemble_face_term2(fe_v_subface_x.get_present_fe_values (),
- fe_v_face_neighbor_x.get_present_fe_values (),
- u_v_matrix,
- un_v_matrix,
- u_vn_matrix,
- un_vn_matrix);
-
- dofs_neighbor.resize (dofs_on_neighbor_child);
- neighbor_child->get_dof_indices (dofs_neighbor);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_on_neighbor_child; ++j)
- system_matrix.add(dofs[i], dofs_neighbor[j],
- un_v_matrix(i,j));
-
- for (unsigned int i=0; i<dofs_on_neighbor_child; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- system_matrix.add(dofs_neighbor[i], dofs[j],
- u_vn_matrix(i,j));
-
- for (unsigned int i=0; i<dofs_on_neighbor_child; ++i)
- for (unsigned int j=0; j<dofs_on_neighbor_child; ++j)
- system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
- un_vn_matrix(i,j));
- }
- }
- else
- {
- // Case 3, with the
- // additional rule
- // a)
- if (neighbor->level() == cell->level() &&
- neighbor->index() > cell->index())
- {
- const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
-
- const unsigned int quadrature_index =
- std::max (neighbor->active_fe_index (),
- cell->active_fe_index ());
-
- un_v_matrix = 0;
- u_vn_matrix = 0;
- un_vn_matrix = 0;
-
- fe_v_face_x.reinit (cell, face_no, quadrature_index);
- fe_v_face_neighbor_x.reinit (neighbor, neighbor2, quadrature_index);
-
- dg.assemble_face_term2(fe_v_face_x.get_present_fe_values (),
- fe_v_face_neighbor_x.get_present_fe_values (),
- u_v_matrix,
- un_v_matrix,
- u_vn_matrix,
- un_vn_matrix);
-
- dofs_neighbor.resize (dofs_on_neighbor);
- neighbor->get_dof_indices (dofs_neighbor);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_on_neighbor; ++j)
- system_matrix.add(dofs[i], dofs_neighbor[j],
- un_v_matrix(i,j));
-
- for (unsigned int i=0; i<dofs_on_neighbor; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- system_matrix.add(dofs_neighbor[i], dofs[j],
- u_vn_matrix(i,j));
-
- for (unsigned int i=0; i<dofs_on_neighbor; ++i)
- for (unsigned int j=0; j<dofs_on_neighbor; ++j)
- system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
- un_vn_matrix(i,j));
-
- }
-
- // Due to rule b)
- // we do not need
- // to consider case
- // 4.
- }
- }
- }
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- system_matrix.add(dofs[i], dofs[j], u_v_matrix(i,j));
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- right_hand_side(dofs[i]) += cell_vector(i);
- }
}
- // @sect3{All the rest}
- //
- // First, we have to solve the
- // discrete system. Since we solve a
- // transport equation, the matrix is
- // nonsymmetric. Hence, we use a
- // GMRES solver.
+
+
+ // @sect3{Equation data}
+
+ // As this is one more academic
+ // example, we'd like to compare
+ // exact and computed solution
+ // against each other. For this, we
+ // need to declare function classes
+ // representing the exact solution
+ // (for comparison and for the
+ // Dirichlet boundary values), as
+ // well as a class that denotes the
+ // right hand side of the equation
+ // (this is simply the Laplace
+ // operator applied to the exact
+ // solution we'd like to recover).
//
- // For a preconditioner, we use the
- // ILU method. Since we already
- // sorted the degrees of freedom in
- // downwind direction, this should be
- // quite efficient. Actually, a block
- // Gauss-Seidel method would be an
- // exact solver, but it has not been
- // implemented for variable block
- // sizes yet.
+ // For this example, let us choose as
+ // exact solution the function
+ // $u(x,y)=exp(x+sin(10y+5x^2))$. In more
+ // than two dimensions, simply repeat
+ // the sine-factor with <code>y</code>
+ // replaced by <code>z</code> and so on. Given
+ // this, the following two classes
+ // are probably straightforward from
+ // the previous examples.
//
+ // As in previous examples, the C++
+ // language forces us to declare and
+ // define a constructor to the
+ // following classes even though they
+ // are empty. This is due to the fact
+ // that the base class has no default
+ // constructor (i.e. one without
+ // arguments), even though it has a
+ // constructor which has default
+ // values for all arguments.
template <int dim>
-void DGMethod<dim>::solve (Vector<double> &solution)
+class Solution : public Function<dim>
{
- SolverControl solver_control (10000, 1e-12, false, true);
- SolverGMRES<Vector<double> > solver (solver_control);
- // Initialize the ILU
- // preconditioner. We decide for
- // two additional off diagonals in
- // order to enhance its
- // performance.
- SparseILU<double>::AdditionalData data(0., 2);
- SparseILU<double> preconditioner;
- preconditioner.initialize (system_matrix, data);
- // Then solve the system:
- solver.solve (system_matrix, solution, right_hand_side,
- preconditioner);
-}
+ public:
+ Solution () : Function<dim> () {};
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const;
+};
- // 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 simpliest
- // 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|$. Futhermore 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 in $H^1$
- // (to be more precise, in $H^1_\beta$)
- // only.
template <int dim>
-void DGMethod<dim>::refine_grid ()
+double
+Solution<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
{
- // 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 (dof_handler,
- solution2,
- gradient_indicator);
-
- // and they are cell-wise scaled by
- // the factor $h^{1+d/2}$
- typename hp::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);
-
- // Simple heuristic hp-Refinement.
- // As the indicator marks nonsmooth
- // regions, p-refine all non marked
- // regions, while the marked
- // regions clearly deserve an
- // h-refinement.
- cell = dof_handler.begin_active ();
- for (; cell!=endc; ++cell)
- if (!cell->refine_flag_set ()
- &&
- (cell->active_fe_index() < fe_collection.size()-1))
- cell->set_active_fe_index (cell->active_fe_index () + 1);
-
- triangulation.execute_coarsening_and_refinement ();
+ double q = p(0);
+ for (unsigned int i=1; i<dim; ++i)
+ q += std::sin(10*p(i)+5*p(0)*p(0));
+ const double exponential = std::exp(q);
+ return exponential;
}
- // 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 DGMethod<dim>::output_results (const unsigned int cycle) const
+class RightHandSide : public Function<dim>
{
- // 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);
+ public:
+ RightHandSide () : Function<dim> () {};
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const;
+};
- Vector<double> active_fe_indices (triangulation.n_active_cells());
- {
- unsigned int index = 0;
- for (typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active();
- cell != dof_handler.end(); ++cell, ++index)
- active_fe_indices(index) = cell->active_fe_index ();
- }
-
-
- // Output of the solution in
- // gnuplot format.
- filename = "sol-";
- filename += ('0' + cycle);
- Assert (cycle < 10, ExcInternalError());
-
- // filename += ".gnuplot";
- filename += ".gmv";
- deallog << "Writing solution to <" << filename << ">..."
- << std::endl << std::endl;
- std::ofstream gnuplot_output (filename.c_str());
-
- DataOut<dim, hp::DoFHandler<dim> > data_out;
- data_out.attach_dof_handler (dof_handler);
- data_out.add_data_vector (solution2, "u");
- data_out.add_data_vector (active_fe_indices, "fe_index");
- data_out.build_patches (4);
-
-// data_out.write_gnuplot(gnuplot_output);
- data_out.write_gmv(gnuplot_output);
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim> &/*p*/,
+ const unsigned int /*component*/) const
+{
+ return 1.;
}
- // The following <code>run</code> function is
- // similar to previous examples. The
- // only difference is that the
- // problem is assembled and solved
- // twice on each refinement step;
- // first by <code>assemble_system1</code> that
- // implements the first version and
- // then by <code>assemble_system2</code> that
- // implements the second version of
- // writing the DG
- // discretization. Furthermore the
- // time needed by each of the two
- // assembling routines is measured.
+
+ // @sect3{The driver routines}
+
+ // What is now missing are only the
+ // functions that actually select the
+ // various options, and run the
+ // simulation on successively finer
+ // grids to monitor the progress as
+ // the mesh is refined.
+ //
+ // This we do in the following
+ // function: it takes a solver
+ // object, and a list of
+ // postprocessing (evaluation)
+ // objects, and runs them with
+ // intermittent mesh refinement:
template <int dim>
-void DGMethod<dim>::run ()
+void
+run_simulation (LaplaceSolver::Base<dim> &solver,
+ const std::list<Evaluation::EvaluationBase<dim> *> &postprocessor_list)
{
- for (unsigned int cycle=0; cycle<7; ++cycle)
+ // We will give an indicator of the
+ // step we are presently computing,
+ // in order to keep the user
+ // informed that something is still
+ // happening, and that the program
+ // is not in an endless loop. This
+ // is the head of this status line:
+ std::cout << "Refinement cycle: ";
+
+ // Then start a loop which only
+ // terminates once the number of
+ // degrees of freedom is larger
+ // than 20,000 (you may of course
+ // change this limit, if you need
+ // more -- or less -- accuracy from
+ // your program).
+ for (unsigned int step=0; true; ++step)
{
- deallog << "Cycle " << cycle << ':' << std::endl;
-
- if (cycle == 0)
+ // Then give the <code>alive</code>
+ // indication for this
+ // iteration. Note that the
+ // <code>std::flush</code> is needed to
+ // have the text actually
+ // appear on the screen, rather
+ // than only in some buffer
+ // that is only flushed the
+ // next time we issue an
+ // end-line.
+ std::cout << step << " " << std::flush;
+
+ // Now solve the problem on the
+ // present grid, and run the
+ // evaluators on it. The long
+ // type name of iterators into
+ // the list is a little
+ // annoying, but could be
+ // shortened by a typedef, if
+ // so desired.
+ solver.solve_problem ();
+
+ for (typename std::list<Evaluation::EvaluationBase<dim> *>::const_iterator
+ i = postprocessor_list.begin();
+ i != postprocessor_list.end(); ++i)
{
- GridGenerator::hyper_cube (triangulation);
+ (*i)->set_refinement_cycle (step);
+ solver.postprocess (**i);
+ };
- triangulation.refine_global (1);
- }
- else
- refine_grid ();
-
- deallog << " Number of active cells: "
- << triangulation.n_active_cells()
- << std::endl;
+ // Now check whether more
+ // iterations are required, or
+ // whether the loop shall be
+ // ended:
+ if (solver.n_dofs() < 20000)
+ solver.refine_grid ();
+ else
+ break;
+ };
- setup_system ();
+ // Finally end the line in which we
+ // displayed status reports:
+ std::cout << std::endl;
+}
- deallog << " Number of degrees of freedom: "
- << dof_handler.n_dofs()
- << std::endl;
- // The constructor of the Timer
- // class automatically starts
- // the time measurement.
- Timer assemble_timer;
- // First assembling routine.
- assemble_system1 ();
- // The operator () accesses the
- // current time without
- // disturbing the time
- // measurement.
- deallog << "Time of assemble_system1: "
- << assemble_timer()
- << std::endl;
- solve (solution1);
-
-
- // As preparation for the
- // second assembling routine we
- // reinit the system matrix, the
- // right hand side vector and
- // the Timer object.
- system_matrix = 0;
- right_hand_side = 0;
- assemble_timer.reset();
-
- // We start the Timer,
- assemble_timer.start();
- // call the second assembling routine
- assemble_system2 ();
- // and access the current time.
- deallog << "Time of assemble_system2: "
- << assemble_timer()
- << std::endl;
- solve (solution2);
-
- // To make sure that both
- // versions of the DG method
- // yield the same
- // discretization and hence the
- // same solution we check the
- // two solutions for equality.
- solution1-=solution2;
-
- const double difference=solution1.linfty_norm();
- if (difference>1e-12)
- deallog << "solution1 and solution2 differ!!" << std::endl;
- else
- deallog << "solution1 and solution2 coincide." << std::endl;
-
- // Finally we perform the
- // output.
- output_results (cycle);
+void
+create_coarse_grid (Triangulation<2> &coarse_grid)
+{
+ const unsigned int dim = 2;
+ static const Point<2> vertices_1[]
+ = { Point<2> (-1., -1.),
+ Point<2> (-1./2, -1.),
+ Point<2> (0., -1.),
+ Point<2> (+1./2, -1.),
+ Point<2> (+1, -1.),
+
+ Point<2> (-1., -1./2.),
+ Point<2> (-1./2, -1./2.),
+ Point<2> (0., -1./2.),
+ Point<2> (+1./2, -1./2.),
+ Point<2> (+1, -1./2.),
+
+ Point<2> (-1., 0.),
+ Point<2> (-1./2, 0.),
+ Point<2> (+1./2, 0.),
+ Point<2> (+1, 0.),
+
+ Point<2> (-1., 1./2.),
+ Point<2> (-1./2, 1./2.),
+ Point<2> (0., 1./2.),
+ Point<2> (+1./2, 1./2.),
+ Point<2> (+1, 1./2.),
+
+ Point<2> (-1., 1.),
+ Point<2> (-1./2, 1.),
+ Point<2> (0., 1.),
+ Point<2> (+1./2, 1.),
+ Point<2> (+1, 1.) };
+ const unsigned int
+ n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]);
+ const std::vector<Point<dim> > vertices (&vertices_1[0],
+ &vertices_1[n_vertices]);
+ static const int cell_vertices[][GeometryInfo<dim>::vertices_per_cell]
+ = {{0, 1, 5, 6},
+ {1, 2, 6, 7},
+ {2, 3, 7, 8},
+ {3, 4, 8, 9},
+ {5, 6, 10, 11},
+ {8, 9, 12, 13},
+ {10, 11, 14, 15},
+ {12, 13, 17, 18},
+ {14, 15, 19, 20},
+ {15, 16, 20, 21},
+ {16, 17, 21, 22},
+ {17, 18, 22, 23}};
+ const unsigned int
+ n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]);
+
+ std::vector<CellData<dim> > cells (n_cells, CellData<dim>());
+ for (unsigned int i=0; i<n_cells; ++i)
+ {
+ for (unsigned int j=0;
+ j<GeometryInfo<dim>::vertices_per_cell;
+ ++j)
+ cells[i].vertices[j] = cell_vertices[i][j];
+ cells[i].material_id = 0;
}
+
+ coarse_grid.create_triangulation (vertices,
+ cells,
+ SubCellData());
+ coarse_grid.refine_global (1);
+}
+
+
+ // The final function is one which
+ // takes the name of a solver
+ // (presently "kelly" and "global"
+ // are allowed), creates a solver
+ // object out of it using a coarse
+ // grid (in this case the ubiquitous
+ // unit square) and a finite element
+ // object (here the likewise
+ // ubiquitous bilinear one), and uses
+ // that solver to ask for the
+ // solution of the problem on a
+ // sequence of successively refined
+ // grids.
+ //
+ // The function also sets up two of
+ // evaluation functions, one
+ // evaluating the solution at the
+ // point (0.5,0.5), the other writing
+ // out the solution to a file.
+template <int dim>
+void solve_problem ()
+{
+ Triangulation<dim> triangulation;
+ create_coarse_grid (triangulation);
+
+ const FE_Q<dim> fe(1);
+ const QGauss<dim> quadrature(4);
+ const RightHandSide<dim> rhs_function;
+ const ZeroFunction<dim> boundary_values;
+
+ // Create a solver object of the
+ // kind indicated by the argument
+ // to this function. If the name is
+ // not recognized, throw an
+ // exception!
+ LaplaceSolver::RefinementKelly<dim> solver (triangulation, fe,
+ quadrature,
+ rhs_function,
+ boundary_values);
+
+ // Next create a table object in
+ // which the values of the
+ // numerical solution at the point
+ // (0.5,0.5) will be stored, and
+ // create a respective evaluation
+ // object:
+ TableHandler results_table;
+ Evaluation::PointValueEvaluation<dim>
+ postprocessor1 (Point<dim>(0.5,0.5), results_table);
+
+ // Also generate an evaluator which
+ // writes out the solution:
+ Evaluation::SolutionOutput<dim>
+ postprocessor2 (std::string("solution"),
+ DataOut<dim>::vtk);
+
+ // Take these two evaluation
+ // objects and put them in a
+ // list...
+ std::list<Evaluation::EvaluationBase<dim> *> postprocessor_list;
+ postprocessor_list.push_back (&postprocessor1);
+ postprocessor_list.push_back (&postprocessor2);
+
+ // ... which we can then pass on to
+ // the function that actually runs
+ // the simulation on successively
+ // refined grids:
+ run_simulation (solver, postprocessor_list);
+
+ // When this all is done, write out
+ // the results of the point
+ // evaluations, and finally delete
+ // the solver object:
+ results_table.write_text (std::cout);
+
+ // And one blank line after all
+ // results:
+ std::cout << std::endl;
}
- // The following <code>main</code> function is
- // similar to previous examples and
- // need not to be commented on.
+
+
+ // There is not much to say about the
+ // main function. It follows the same
+ // pattern as in all previous
+ // examples, with attempts to catch
+ // thrown exceptions, and displaying
+ // as much information as possible if
+ // we should get some. The rest is
+ // self-explanatory.
int main ()
{
try
{
- DGMethod<2> dgmethod;
- dgmethod.run ();
+ deallog.depth_console (0);
+
+ solve_problem<2> ();
}
catch (std::exception &exc)
{
<< std::endl;
return 1;
};
-
+
return 0;
}
-
-