From: bangerth Date: Wed, 13 Dec 2006 05:44:44 +0000 (+0000) Subject: Replace the hp adaptation of step-12 by a mixture between step-13 and step-14. It... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=af0a0b06db43f223b9a49730e6ab4ff64d8783d4;p=dealii-svn.git Replace the hp adaptation of step-12 by a mixture between step-13 and step-14. It still has to be hp-ified, but has more promise and is less of a crutch of a code. git-svn-id: https://svn.dealii.org/trunk@14229 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc index 37cc153d32..8cca3c9577 100644 --- a/deal.II/examples/step-27/step-27.cc +++ b/deal.II/examples/step-27/step-27.cc @@ -1,1786 +1,2034 @@ /* $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 + // base -- lac -- grid -- + // dofs -- fe -- numerics + // (as each of these categories + // roughly builds upon previous + // ones), then C++ standard headers: #include #include +#include +#include +#include #include +#include #include -#include -#include +#include #include -#include #include #include -#include -#include #include #include -#include - // 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 - // 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 +#include +#include #include #include -#include +#include +#include +#include +#include #include +#include -#include -#include - - // We are going to use gradients as - // refinement indicator. -#include - // Finally we do some time comparison - // using the Timer class. -#include - - // And this again is C++: + // Now for the C++ standard headers: #include #include - - // 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 -#include - - // 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 - - // 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 - +#include +#include // 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 + // DoFHandler 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, RHS - // and BoundaryValues, are - // derived from the Function - // class. Only the value_list - // function are implemented because - // only lists of function values are - // computed rather than single - // values. -template -class RHS: public Function + // From an abstract point of view, we + // declare a pure base class + // that provides an evaluation + // operator operator() which will + // do the evaluation of the solution + // (whatever derived classes might + // consider an evaluation). 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 operator() a + // functor 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 > &points, - std::vector &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 operator() + // taking a DoFHandler 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 + class EvaluationBase + { + public: + virtual ~EvaluationBase (); -template -class BoundaryValues: public Function -{ - public: - virtual void value_list (const std::vector > &points, - std::vector &values, - const unsigned int component=0) const; -}; + void set_refinement_cycle (const unsigned int refinement_cycle); + + virtual void operator () (const DoFHandler &dof_handler, + const Vector &solution) const = 0; + protected: + unsigned int refinement_cycle; + }; + + + // After the declaration has been + // discussed above, the + // implementation is rather + // straightforward: + template + EvaluationBase::~EvaluationBase () + {} + + + template + void + EvaluationBase::set_refinement_cycle (const unsigned int step) + { + refinement_cycle = step; + } - // The class Beta represents the - // vector valued flow field of the - // linear transport equation and is - // not derived from the Function - // class as we prefer to get function - // values of type Point rather - // than of type - // Vector@. This, because - // there exist scalar products - // between Point and Point as - // well as between Point and - // Tensor, 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 -class Beta -{ - public: - Beta () {}; - void value_list (const std::vector > &points, - std::vector > &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 + // operator(). 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 DeclExceptionN + // 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 + class PointValueEvaluation : public EvaluationBase + { + public: + PointValueEvaluation (const Point &evaluation_point, + TableHandler &results_table); + + virtual void operator () (const DoFHandler &dof_handler, + const Vector &solution) const; + + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + private: + const Point evaluation_point; + TableHandler &results_table; + }; + + + // As for the definition, the + // constructor is trivial, just + // taking data and storing it in + // object-local ones: + template + PointValueEvaluation:: + PointValueEvaluation (const Point &evaluation_point, + TableHandler &results_table) + : + evaluation_point (evaluation_point), + results_table (results_table) + {} + - // The implementation of the - // value_list functions of these - // classes are rather simple. For - // simplicity the right hand side is - // set to be zero but will be - // assembled anyway. -template -void RHS::value_list(const std::vector > &points, - std::vector &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 + void + PointValueEvaluation:: + operator () (const DoFHandler &dof_handler, + const Vector &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::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::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 + // Assert + // (dof_handler.get_fe().dofs_per_vertex + // @> 0, + // ExcNotImplemented()), + // 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 + // cell-@>vertex_dof_index(vertex,0) + // 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 AssertThrow 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 throw + // 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 catch 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 + // Assert macro in other + // example programs as well. It + // differed from the + // AssertThrow 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 AssertThrow + // 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 -void Beta::value_list(const std::vector > &points, - std::vector > &values) const -{ - Assert(values.size()==points.size(), - ExcDimensionMismatch(values.size(),points.size())); - for (unsigned int i=0; i &p=points[i]; - Point &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 evaluation of a + // solution is to output it to a + // file in a graphical + // format. Since in the evaluation + // functions we are given a + // DoFHandler 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 EvaluationBase base + // class, its main interface is the + // operator() + // 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 + // DataOutInterface class + // (which is a base class of + // DataOut through which we + // will access its fields) provides + // an enumeration field + // OutputFormat, 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 ucd, + // gnuplot, povray, + // eps, gmv, tecplot, + // tecplot_binary, dx, and + // vtk, 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 + // write function, which then + // branches to the + // write_gnuplot, + // write_ucd, 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 + class SolutionOutput : public EvaluationBase + { + public: + SolutionOutput (const std::string &output_name_base, + const typename DataOut::OutputFormat output_format); + + virtual void operator () (const DoFHandler &dof_handler, + const Vector &solution) const; + private: + const std::string output_name_base; + const typename DataOut::OutputFormat output_format; + }; + + + template + SolutionOutput:: + SolutionOutput (const std::string &output_name_base, + const typename DataOut::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 + // DataOut::default_suffix + // 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 + // DataOut::write 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 + // this-@> 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 two-stage + // name lookup there). + template + void + SolutionOutput::operator () (const DoFHandler &dof_handler, + const Vector &solution) const + { + DataOut 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 -void BoundaryValues::value_list(const std::vector > &points, - std::vector &values, - const unsigned int) const -{ - Assert(values.size()==points.size(), - ExcDimensionMismatch(values.size(),points.size())); - for (unsigned int i=0; iDGTransportEquation. Its - // member functions were already - // mentioned in the Introduction and - // will be explained - // below. Furthermore it includes - // objects of the previously defined - // Beta, RHS and - // BoundaryValues function - // classes. -template -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& fe_v, - FullMatrix &u_v_matrix, - Vector &cell_vector) const; - - void assemble_boundary_term(const FEFaceValues& fe_v, - FullMatrix &u_v_matrix, - Vector &cell_vector) const; - - void assemble_face_term1(const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &u_v_matrix, - FullMatrix &un_v_matrix) const; - - void assemble_face_term2(const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &u_v_matrix, - FullMatrix &un_v_matrix, - FullMatrix &u_vn_matrix, - FullMatrix &un_vn_matrix) const; - private: - const Beta beta_function; - const RHS rhs_function; - const BoundaryValues 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 Object, 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 SmartPointer 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 + // refine_grid function. + // + // Finally, we have a function + // n_dofs 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 + class Base + { + public: + Base (Triangulation &coarse_grid); + virtual ~Base (); + + virtual void solve_problem () = 0; + virtual void postprocess (const Evaluation::EvaluationBase &postprocessor) const = 0; + virtual void refine_grid () = 0; + virtual unsigned int n_dofs () const = 0; + + protected: + const SmartPointer > triangulation; + }; -template -DGTransportEquation::DGTransportEquation () - : - beta_function (), - rhs_function (), - boundary_function () -{} + // The implementation of the only + // two non-abstract functions is + // then rather boring: + template + Base::Base (Triangulation &coarse_grid) + : + triangulation (&coarse_grid) + {} - // @sect4{Function: assemble_cell_term} - // - // The assemble_cell_term - // function assembles the cell terms - // of the discretization. - // u_v_matrix is a cell matrix, - // i.e. for a DG method of degree 1, - // it is of size 4 times 4, and - // cell_vector is of size 4. - // When this function is invoked, - // fe_v is already reinit'ed with the - // current cell before and includes - // all shape values needed. -template -void DGTransportEquation::assemble_cell_term( - const FEValues &fe_v, - FullMatrix &u_v_matrix, - Vector &cell_vector) const -{ - // First we ask fe_v for the - // shape gradients, shape values and - // quadrature weights, - const std::vector &JxW = fe_v.get_JxW_values (); - - // Then the flow field beta and the - // rhs_function are evaluated at - // the quadrature points, - std::vector > beta (fe_v.n_quadrature_points); - std::vector 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 + Base::~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; pointsolve_problem and + // postprocess functions + // declared in the base class. It + // does not, however, implement the + // refine_grid method, as mesh + // refinement will be implemented + // in a number of derived classes. + // + // It also declares a new abstract + // virtual function, + // assemble_rhs, 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 + // DoFHandler object + // later. Since finite elements and + // quadrature formula should match, + // it is also passed a quadrature + // object. + // + // The solve_problem sets up + // the data structures for the + // actual solution, calls the + // functions to assemble the linear + // system, and solves it. + // + // The postprocess function + // finally takes an evaluation + // object and applies it to the + // computed solution. + // + // The n_dofs function finally + // implements the pure virtual + // function of the base class. + template + class Solver : public virtual Base + { + public: + Solver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &boundary_values); + virtual + ~Solver (); + + virtual + void + solve_problem (); + + virtual + void + postprocess (const Evaluation::EvaluationBase &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 > fe; + const SmartPointer > quadrature; + DoFHandler dof_handler; + Vector solution; + const SmartPointer > 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 &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 &dof_handler); + + void solve (Vector &solution) const; - cell_vector(i) += rhs[point] * fe_v.shape_value(i,point) * JxW[point]; - } -} + ConstraintMatrix hanging_node_constraints; + SparsityPattern sparsity_pattern; + SparseMatrix matrix; + Vector 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::active_cell_iterator &begin_cell, + const typename DoFHandler::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 + // DoFHandler 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 + // solve_problem function). + template + Solver::Solver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &boundary_values) + : + Base (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 + Solver::~Solver () + { + dof_handler.clear (); + } - // @sect4{Function: assemble_boundary_term} - // - // The assemble_boundary_term - // function assembles the face terms - // at boundary faces. When this - // function is invoked, fe_v is - // already reinit'ed with the current - // cell and current face. Hence it - // provides the shape values on that - // boundary face. -template -void DGTransportEquation::assemble_boundary_term( - const FEFaceValues& fe_v, - FullMatrix &u_v_matrix, - Vector &cell_vector) const -{ - // Again, as in the previous - // function, we ask the FEValues - // object for the shape values and - // the quadrature weights - const std::vector &JxW = fe_v.get_JxW_values (); - // but here also for the normals. - const std::vector > &normals = fe_v.get_normal_vectors (); - - // We evaluate the flow field - // and the boundary values at the - // quadrature points. - std::vector > beta (fe_v.n_quadrature_points); - std::vector 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; point0) - for (unsigned int i=0; i + void + Solver::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 assemble_face_term1 - // 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, - // fe_v and fe_v_neighbor 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 - // u_v_matrix this function has - // got a new argument - // un_v_matrix, 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 - // un is the short notation for - // u_neighbor and represents - // $\hat u_h$. -template -void DGTransportEquation::assemble_face_term1( - const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &u_v_matrix, - FullMatrix &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 &JxW = fe_v.get_JxW_values (); - const std::vector > &normals = fe_v.get_normal_vectors (); - - // and we evaluate the flow field - // at the quadrature points. - std::vector > 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; point0) - for (unsigned int i=0; iassemble_face_term2 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 - // u_vn_matrix and - // un_vn_matrix that will store - // contributions due to terms - // involving u and vn as well as un - // and vn. -template -void DGTransportEquation::assemble_face_term2( - const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &u_v_matrix, - FullMatrix &un_v_matrix, - FullMatrix &u_vn_matrix, - FullMatrix &un_vn_matrix) const -{ - // the first few lines are the same - const std::vector &JxW = fe_v.get_JxW_values (); - const std::vector > &normals = fe_v.get_normal_vectors (); + assemble_linear_system (linear_system); + linear_system.solve (solution); + } - std::vector > beta (fe_v.n_quadrature_points); - - beta_function.value_list (fe_v.get_quadrature_points(), beta); - for (unsigned int point=0; point0) - { - // This terms we've already seen. - for (unsigned int i=0; ipostprocess 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 + void + Solver:: + postprocess (const Evaluation::EvaluationBase &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 DGMethod 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 -class DGMethod -{ - public: - DGMethod (); - ~DGMethod (); + // The n_dofs function should + // be self-explanatory: + template + unsigned int + Solver::n_dofs () const + { + return dof_handler.n_dofs(); + } + - void run (); - - private: - void setup_system (); - void assemble_system1 (); - void assemble_system2 (); - void solve (Vector &solution); - void refine_grid (); - void output_results (const unsigned int cycle) const; - - Triangulation 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 + void + Solver::assemble_linear_system (LinearSystem &linear_system) + { + // First define a convenience + // abbreviation for these lengthy + // iterator names... + typedef + typename DoFHandler::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 > + thread_ranges + = Threads::split_range (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::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 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 fe_collection; - - // As already mentioned, the - // standard DoFHandler has to be - // replaced by a hp::DoFHandler. - hp::DoFHandler dof_handler; - - SparsityPattern sparsity; - SparseMatrix 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 quadratures; - hp::QCollection face_quadratures; - // And there are two solution - // vectors, that store the - // solutions to the problems - // corresponding to the two - // different assembling routines - // assemble_system1 and - // assemble_system2; - Vector solution1; - Vector solution2; - Vector 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 + void + Solver::assemble_matrix (LinearSystem &linear_system, + const typename DoFHandler::active_cell_iterator &begin_cell, + const typename DoFHandler::active_cell_iterator &end_cell, + Threads::ThreadMutex &mutex) const + { + FEValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + for (typename DoFHandler::active_cell_iterator cell=begin_cell; + cell!=end_cell; ++cell) + { + cell_matrix = 0; + + fe_values.reinit (cell); + + for (unsigned int q_point=0; q_pointget_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 acquire and + // release 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 + // return, 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 + // scoped lock 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 + // return 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; ilock 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 Threads::ThreadGroup + // 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 + // DoFTools::make_hanging_node_constraints + // 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 & 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 mhnc_p (short for + // pointer to + // make_hanging_node_constraints) + // with the right type, and using + // this pointer instead. + template + Solver::LinearSystem:: + LinearSystem (const DoFHandler &dof_handler) + { + hanging_node_constraints.clear (); + + void (*mhnc_p) (const DoFHandler &, + ConstraintMatrix &) + = &DoFTools::make_hanging_node_constraints; - // Finally this class includes an - // object of the - // DGTransportEquations class - // described above. - const DGTransportEquation 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 + // hanging_node_constraints + // 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 -DGMethod::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 (i)); - quadratures.push_back (QGauss (i+2)); - face_quadratures.push_back (QGauss (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 + void + Solver::LinearSystem::solve (Vector &solution) const + { + SolverControl solver_control (1000, 1e-12); + SolverCG<> cg (solver_control); -template -DGMethod::~DGMethod () -{ - dof_handler.clear (); -} + PreconditionSSOR<> preconditioner; + preconditioner.initialize(matrix, 1.2); + cg.solve (matrix, solution, rhs, preconditioner); -template -void DGMethod::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 sorting_direction; - for (unsigned int d=0;dassemble_system1 function that - // implements the DG discretization - // in its first version. This - // function repeatedly calls the - // assemble_cell_term, - // assemble_boundary_term and - // assemble_face_term1 functions - // of the DGTransportEquation - // object. The - // assemble_boundary_term covers - // the first case mentioned in the - // introduction. - // - // 1. face is at boundary - // - // This function takes a - // FEFaceValues object as - // argument. In contrast to that - // assemble_face_term1 - // takes two FEFaceValuesBase - // 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 FEFaceValues - // or of class FESubfaceValues - // (both derived from - // FEFaceValuesBase) according to - // the remaining cases mentioned - // in the introduction: - // - // 2. neighboring cell is finer - // (current cell: FESubfaceValues, - // neighboring cell: FEFaceValues); - // - // 3. neighboring cell is of the same - // refinement level (both, current - // and neighboring cell: - // FEFaceValues); - // - // 4. neighboring cell is coarser - // (current cell: FEFaceValues, - // neighboring cell: - // FESubfaceValues). - // - // 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 -void DGMethod::assemble_system1 () -{ - // First we create the - // UpdateFlags for the - // FEValues and the - // FEFaceValues 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 FEValues - // object. Here, we use the default - // MappingQ1. different mapping - // create a MappingCollection first - // and call the respective - // hp::FEValues constructor. - hp::FEValues fe_v_x (fe_collection, quadratures, update_flags); - - // Similarly we create the - // FEFaceValues and - // FESubfaceValues 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 fe_v_face_x ( - fe_collection, face_quadratures, face_update_flags); - hp::FESubfaceValues fe_v_subface_x ( - fe_collection, face_quadratures, face_update_flags); - hp::FEFaceValues fe_v_face_neighbor_x ( - fe_collection, face_quadratures, neighbor_face_update_flags); - hp::FESubfaceValues 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 v (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 u_v_matrix (max_dofs_per_cell, max_dofs_per_cell); - FullMatrix un_v_matrix (max_dofs_per_cell, max_dofs_per_cell); - Vector cell_vector (max_dofs_per_cell); - - // Furthermore we need some cell - // iterators. - typename hp::DoFHandler::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 dofs (dofs_per_cell); - std::vector dofs_neighbor; - - // In the - // assemble_face_term1 - // function contributions to - // the cell matrices and the - // cell vector are only - // ADDED. Therefore on each - // cell we need to reset the - // u_v_matrix and - // cell_vector to zero, - // before assembling the cell terms. - u_v_matrix = 0; - cell_vector = 0; - - // Now we reinit the FEValues - // object for the current cell - fe_v_x.reinit (cell); - - // and call the function - // that assembles the cell - // terms. The first argument is - // the FEValues 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::faces_per_cell; ++face_no) - { - // First we set the face - // iterator - typename hp::DoFHandler::face_iterator face=cell->face(face_no); - - // and clear the - // un_v_matrix 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 - // FEFaceValues - // 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::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 dim==1. - 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 - // cell. - const unsigned int neighbor2= - cell->neighbor_of_neighbor(face_no); - - - // We loop over - // subfaces - for (unsigned int subface_no=0; - subface_non_children(); ++subface_no) - { - // and set the - // cell - // iterator - // neighbor_child - // to the cell - // placed - // `behind' the - // current - // subface. - typename hp::DoFHandler::active_cell_iterator neighbor_child= - neighbor->child(GeometryInfo:: - 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 - // un_v_matrix - // on each - // subface - // because on - // each subface - // the un - // belong to - // different - // neighboring - // cells. - un_v_matrix = 0; - - // As already - // mentioned - // above for - // the current - // case (case - // 2) we employ - // the - // FESubfaceValues - // 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 - // un_v_matrix - // to the - // system_matrix - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) - for (unsigned int k=0; kget_fe().dofs_per_cell; ++k) - system_matrix.add(dofs[i], dofs_neighbor[k], - un_v_matrix(i,k)); - } - // End of if - // (face-@>has_children()) - } - 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 - // FEFaceValues - // 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 if - // (neighbor-@>level() - // == - // cell-@>level()) - } - 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 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::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 - // FEFaceValues - // 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 - // neighbor_child - // cell, - dofs_neighbor.resize (neighbor->get_fe().dofs_per_cell); - neighbor->get_dof_indices (dofs_neighbor); - - // and distribute the - // un_v_matrix. - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) - for (unsigned int k=0; kget_fe().dofs_per_cell; ++k) - system_matrix.add(dofs[i], dofs_neighbor[k], - un_v_matrix(i,k)); - } - // End of face not at boundary: - } - // End of loop over all faces: - } - - // Finally we distribute the - // u_v_matrix - for (unsigned int i=0; iSmartPointer, 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 assemble_rhs + // method that does what its name + // suggests. + template + class PrimalSolver : public Solver + { + public: + PrimalSolver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &rhs_function, + const Function &boundary_values); + protected: + const SmartPointer > rhs_function; + virtual void assemble_rhs (Vector &rhs) const; + }; + + + // The constructor of this class + // basically does what it is + // announced to do above... + template + PrimalSolver:: + PrimalSolver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &rhs_function, + const Function &boundary_values) + : + Base (triangulation), + Solver (triangulation, fe, + quadrature, boundary_values), + rhs_function (&rhs_function) + {} + + + + // ... as does the assemble_rhs + // function. Since this is + // explained in several of the + // previous example programs, we + // leave it at that. + template + void + PrimalSolver:: + assemble_rhs (Vector &rhs) const + { + FEValues 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 cell_rhs (dofs_per_cell); + std::vector rhs_values (n_q_points); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::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; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; iassemble_system2 function that - // implements the DG discretization - // in its second version. This - // function is very similar to the - // assemble_system1 - // function. Therefore, here we only - // discuss the differences between - // the two functions. This function - // repeatedly calls the - // assemble_face_term2 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 + class RefinementKelly : public PrimalSolver + { + public: + RefinementKelly (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &rhs_function, + const Function &boundary_values); -template -void DGMethod::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 - // fe_v_face_neighbor as case 4 - // does not occur. - hp::FEValues fe_v_x ( - fe_collection, quadratures, update_flags); - hp::FEFaceValues fe_v_face_x ( - fe_collection, face_quadratures, face_update_flags); - hp::FESubfaceValues fe_v_subface_x ( - fe_collection, face_quadratures, face_update_flags); - hp::FEFaceValues 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 u_v_matrix (max_dofs_per_cell, max_dofs_per_cell); - FullMatrix 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 vn (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 u_vn_matrix (max_dofs_per_cell, max_dofs_per_cell); - FullMatrix un_vn_matrix (max_dofs_per_cell, max_dofs_per_cell); - - Vector cell_vector (max_dofs_per_cell); - - // The following lines are roughly - // the same as in the previous - // function. - typename hp::DoFHandler::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 dofs (dofs_per_cell); - std::vector 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 + RefinementKelly:: + RefinementKelly (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &rhs_function, + const Function &boundary_values) + : + Base (coarse_grid), + PrimalSolver (coarse_grid, fe, quadrature, + rhs_function, boundary_values) + {} + + + + template + void + RefinementKelly::refine_grid () + { + Vector estimated_error_per_cell (this->triangulation->n_active_cells()); + KellyErrorEstimator::estimate (this->dof_handler, + QGauss(3), + typename FunctionMap::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::faces_per_cell; ++face_no) - { - typename hp::DoFHandler::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::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_non_children(); ++subface_no) - { - typename hp::DoFHandler::cell_iterator neighbor_child= - neighbor->child(GeometryInfo::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; ilevel() == 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; iy + // replaced by z 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 -void DGMethod::solve (Vector &solution) +class Solution : public Function { - SolverControl solver_control (10000, 1e-12, false, true); - SolverGMRES > solver (solver_control); - // Initialize the ILU - // preconditioner. We decide for - // two additional off diagonals in - // order to enhance its - // performance. - SparseILU::AdditionalData data(0., 2); - SparseILU preconditioner; - preconditioner.initialize (system_matrix, data); - // Then solve the system: - solver.solve (system_matrix, solution, right_hand_side, - preconditioner); -} + public: + Solution () : Function () {}; + + virtual double value (const Point &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 - // DerivativeApproximation class - // that computes the approximate - // gradients in a way similar to the - // GradientEstimation described - // in Step 9 of this tutorial. In - // fact, the - // DerivativeApproximation class - // was developed following the - // GradientEstimation class of - // Step 9. Relating to the - // discussion in Step 9, here we - // consider $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 -void DGMethod::refine_grid () +double +Solution::value (const Point &p, + const unsigned int /*component*/) const { - // The DerivativeApproximation - // class computes the gradients to - // float precision. This is - // sufficient as they are - // approximate and serve as - // refinement indicators only. - Vector 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::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 -void DGMethod::output_results (const unsigned int cycle) const +class RightHandSide : public Function { - // 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 () {}; + + virtual double value (const Point &p, + const unsigned int component) const; +}; - Vector active_fe_indices (triangulation.n_active_cells()); - { - unsigned int index = 0; - for (typename hp::DoFHandler::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 > 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 +double +RightHandSide::value (const Point &/*p*/, + const unsigned int /*component*/) const +{ + return 1.; } - // The following run function is - // similar to previous examples. The - // only difference is that the - // problem is assembled and solved - // twice on each refinement step; - // first by assemble_system1 that - // implements the first version and - // then by assemble_system2 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 -void DGMethod::run () +void +run_simulation (LaplaceSolver::Base &solver, + const std::list *> &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 alive + // indication for this + // iteration. Note that the + // std::flush 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 *>::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 > vertices (&vertices_1[0], + &vertices_1[n_vertices]); + static const int cell_vertices[][GeometryInfo::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 > cells (n_cells, CellData()); + for (unsigned int i=0; i::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 +void solve_problem () +{ + Triangulation triangulation; + create_coarse_grid (triangulation); + + const FE_Q fe(1); + const QGauss quadrature(4); + const RightHandSide rhs_function; + const ZeroFunction 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 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 + postprocessor1 (Point(0.5,0.5), results_table); + + // Also generate an evaluator which + // writes out the solution: + Evaluation::SolutionOutput + postprocessor2 (std::string("solution"), + DataOut::vtk); + + // Take these two evaluation + // objects and put them in a + // list... + std::list *> 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 main 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) { @@ -1805,8 +2053,6 @@ int main () << std::endl; return 1; }; - + return 0; } - -