From 9bf25f6395390ce4049fdcbeacd7c08dc379cae7 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 18 Mar 2013 23:11:22 +0000 Subject: [PATCH] Add a heat equation solver. Remove the file with surface points that belongs to the previous incarnation of step-26 that was never finished. git-svn-id: https://svn.dealii.org/trunk@28935 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-26/step-26.cc | 754 +++++++--------------------- 1 file changed, 192 insertions(+), 562 deletions(-) diff --git a/deal.II/examples/step-26/step-26.cc b/deal.II/examples/step-26/step-26.cc index 25ebd7f654..136d4aecb2 100644 --- a/deal.II/examples/step-26/step-26.cc +++ b/deal.II/examples/step-26/step-26.cc @@ -1,672 +1,302 @@ -/* Author: Wolfgang Bangerth, University of Heidelberg, 1999 */ +/* Author: Wolfgang Bangerth, Texas A&M University, 2008 */ /* $Id$ */ /* */ -/* Copyright (C) 1999-2007, 2011-2012 by the deal.II authors */ +/* Copyright (C) 2013 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. */ -// @sect3{Include files} - -// The first few (many?) include files have already been used in the previous -// example, so we will not explain their meaning here again. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include -#include -#include #include +#include #include #include +#include #include -#include #include +#include #include -#include -#include +#include #include #include -#include -#include -#include -#include -#include -#include #include #include -// This is new, however: in the previous example we got some unwanted output -// from the linear solvers. If we want to suppress it, we have to include this -// file and add a single line somewhere to the program (see the main() -// function below for that): -#include - -#include -#include - -// The last step is as in all previous programs: namespace Step26 { using namespace dealii; - class PointCloudSurface : public StraightBoundary<3> + template + class HeatEquation { public: - /** - * Constructor. - */ - PointCloudSurface (const std::string &filename); - - /** - * Let the new point be the - * arithmetic mean of the two - * vertices of the line. - * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class for more - * information. - */ - virtual Point<3> - get_new_point_on_line (const Triangulation<3>::line_iterator &line) const; - - /** - * Let the new point be the - * arithmetic mean of the four - * vertices of this quad and the - * four midpoints of the lines, - * which are already created at - * the time of calling this - * function. - * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class for more - * information. - */ - virtual Point<3> - get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const; - - /** - * Gives n=points.size() - * points that splits the - * StraightBoundary line into - * $n+1$ partitions of equal - * lengths. - * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - */ - virtual void - get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line, - std::vector > &points) const; - - /** - * Gives n=points.size()=m*m - * points that splits the - * p{StraightBoundary} quad into - * (m+1)(m+1) subquads of equal - * size. - * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - */ - virtual void - get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, - std::vector > &points) const; - - /** - * A function that, given a point - * p, returns the closest - * point on the surface defined by the - * input file. For the time being, we - * simply return the closest point in the - * point cloud, rather than doing any - * sort of interpolation. - */ - Point<3> closest_point (const Point<3> &p) const; - private: - std::vector > point_list; - }; - - - PointCloudSurface::PointCloudSurface (const std::string &filename) - { - // first read in all the points - { - std::ifstream in (filename.c_str()); - AssertThrow (in, ExcIO()); - - while (in) - { - Point<3> p; - in >> p; - point_list.push_back (p); - } + HeatEquation(); + void run(); - AssertThrow (point_list.size() > 1, ExcIO()); - } + private: + void setup_system(); + void solve_u(); + void output_results() const; - // next fit a linear model through the data cloud to rectify it in a local - // coordinate system - // - // the first step is to move the center of mass of the points to the - // origin - { - const Point<3> c_o_m = std::accumulate (point_list.begin(), - point_list.end(), - Point<3>()) / - point_list.size(); - for (unsigned int i=0; i triangulation; + FE_Q fe; + DoFHandler dof_handler; - // next do a least squares fit to the function ax+by. this leads to the - // following equations: + ConstraintMatrix constraints; - // min f(a,b) = sum_i (zi-a xi - b yi)^2 / 2 - // - // f_a = sum_i (zi - a xi - b yi) xi = 0 f_b = sum_i (zi - a xi - b yi) yi - // = 0 - // - // f_a = (sum_i zi xi) - (sum xi^2) a - (sum xi yi) b = 0 f_a = (sum_i zi - // yi) - (sum xi yi) a - (sum yi^2) b = 0 - { - double A[2][2] = {{0,0},{0,0}}; - double B[2] = {0,0}; + SparsityPattern sparsity_pattern; + SparseMatrix mass_matrix; + SparseMatrix laplace_matrix; + SparseMatrix matrix_u; - for (unsigned int i=0; i solution_u; + Vector old_solution_u; + Vector system_rhs; - B[0] += point_list[i][0] * point_list[i][2]; - B[1] += point_list[i][1] * point_list[i][2]; - } + double time, time_step; + unsigned int timestep_number; + const double theta; + }; - const double det = A[0][0]*A[1][1]-2*A[0][1]; - const double a = (A[1][1] * B[0] - A[0][1] * B[1]) / det; - const double b = (A[0][0] * B[1] - A[0][1] * B[0]) / det; + //------------------------------------- + template + class RightHandSide: public Function + { + public: + RightHandSide() : + Function(), + period (0.2) + {} - // with this information, we can rotate the points so that the - // corresponding least-squares fit would be the x-y plane - const Point<2> gradient_direction - = Point<2>(a,b) / std::sqrt(a*a+b*b); - const Point<2> orthogonal_direction - = Point<2>(-b,a) / std::sqrt(a*a+b*b); + virtual double value(const Point &p, + const unsigned int component = 0) const; - const double stretch_factor = std::sqrt(1.+a*a+b*b); + private: + const double period; + }; - for (unsigned int i=0; i xy (point_list[i][0], - point_list[i][1]); - const double grad_distance = xy * gradient_direction; - const double orth_distance = xy * orthogonal_direction; - - // we then have to stretch the points in the gradient direction. the - // stretch factor is defined above (zero if the original plane was - // already the xy plane, infinity if it was vertical) - const Point<2> new_xy - = (grad_distance * stretch_factor * gradient_direction + - orth_distance * orthogonal_direction); - point_list[i][0] = new_xy[0]; - point_list[i][1] = new_xy[1]; - } - } - } + template + double RightHandSide::value(const Point &p, + const unsigned int component) const + { + Assert (component == 0, ExcInternalError()); + Assert (dim == 2, ExcNotImplemented()); + return std::cos(p[0]*numbers::PI/2) * std::cos(p[1]*numbers::PI/2); - Point<3> - PointCloudSurface::closest_point (const Point<3> &p) const - { - double distance = p.distance (point_list[0]); - Point<3> point = point_list[0]; + const double time = this->get_time(); + const double point_within_period = (time/period - std::floor(time/period)); - for (std::vector >::const_iterator i=point_list.begin(); - i != point_list.end(); ++i) + if ((point_within_period >= 0.0) && (point_within_period <= 0.2)) { - const double d = p.distance (*i); - if (d < distance) - { - distance = d; - point = *i; - } + if ((p[0] > 0.5) && (p[1] > -0.5)) + return 1; + else + return 0; } - - return point; + else if ((point_within_period >= 0.5) && (point_within_period <= 0.7)) + { + if ((p[0] > -0.5) && (p[1] > 0.5)) + return 1; + else + return 0; + } + else + return 0; } - - Point<3> - PointCloudSurface:: - get_new_point_on_line (const Triangulation<3>::line_iterator &line) const + template + class BoundaryValuesU: public Function { - return closest_point (StraightBoundary<3>::get_new_point_on_line (line)); - } - - + public: + BoundaryValuesU() : + Function() + { + } + virtual ~BoundaryValuesU() + { + } + virtual double value(const Point &p, + const unsigned int component = 0) const; + }; - Point<3> - PointCloudSurface:: - get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const + template + double BoundaryValuesU::value(const Point &/*p*/, + const unsigned int component) const { - return closest_point (StraightBoundary<3>::get_new_point_on_quad (quad)); + Assert(component == 0, ExcInternalError()); + return 0; // Zero-Dirichlet Boundary } - - - void - PointCloudSurface:: - get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line, - std::vector > &points) const + template + HeatEquation::HeatEquation() : + fe(1), dof_handler(triangulation), time_step(1. / 1000), theta(0.5) { - StraightBoundary<3>::get_intermediate_points_on_line (line, - points); - for (unsigned int i=0; i::quad_iterator &quad, - std::vector > &points) const + template + void HeatEquation::setup_system() { - StraightBoundary<3>::get_intermediate_points_on_quad (quad, - points); - for (unsigned int i=0; i(3), mass_matrix); + MatrixCreator::create_laplace_matrix(dof_handler, QGauss(3), + laplace_matrix); + solution_u.reinit(dof_handler.n_dofs()); + old_solution_u.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); - // @sect3{The LaplaceProblem class template} + constraints.close(); + } - // This is again the same LaplaceProblem class as in the - // previous example. The only difference is that we have now declared it as - // a class with a template parameter, and the template parameter is of - // course the spatial dimension in which we would like to solve the Laplace - // equation. Of course, several of the member variables depend on this - // dimension as well, in particular the Triangulation class, which has to - // represent quadrilaterals or hexahedra, respectively. Apart from this, - // everything is as before. - template - class LaplaceProblem + template + void HeatEquation::solve_u() { - public: - LaplaceProblem (); - void run (); - - private: - void make_grid_and_dofs (); - void assemble_system (); - void solve (); - void output_results () const; - - Triangulation triangulation; - FE_Q fe; - DoFHandler dof_handler; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; - Vector system_rhs; - }; - - - // @sect3{Right hand side and boundary values} - + SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm()); + SolverCG<> cg(solver_control); + cg.solve(matrix_u, solution_u, system_rhs, PreconditionIdentity()); + std::cout << " u-equation: " << solver_control.last_step() + << " CG iterations." << std::endl; + } - template - class BoundaryValues : public Function + template + void HeatEquation::output_results() const { - public: - BoundaryValues () : Function() {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - }; + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution_u, "U"); + data_out.build_patches(); - template - double BoundaryValues::value (const Point &p, - const unsigned int /*component*/) const - { - return std::max(p[dim-1], -5.); + const std::string filename = "solution-" + + Utilities::int_to_string(timestep_number, 3) + ".vtk"; + std::ofstream output(filename.c_str()); + data_out.write_vtk(output); } - - - // @sect3{Implementation of the LaplaceProblem class} - - // Next for the implementation of the class template that makes use of the - // functions above. As before, we will write everything as templates that - // have a formal parameter dim that we assume unknown at the - // time we define the template functions. Only later, the compiler will find - // a declaration of LaplaceProblem@<2@> (in the - // main function, actually) and compile the entire class with - // dim replaced by 2, a process referred to as `instantiation - // of a template'. When doing so, it will also replace instances of - // RightHandSide@ by RightHandSide@<2@> and - // instantiate the latter class from the class template. - // - // In fact, the compiler will also find a declaration - // LaplaceProblem@<3@> in main(). This will cause - // it to again go back to the general LaplaceProblem@ - // template, replace all occurrences of dim, this time by 3, - // and compile the class a second time. Note that the two instantiations - // LaplaceProblem@<2@> and LaplaceProblem@<3@> are - // completely independent classes; their only common feature is that they - // are both instantiated from the same general template, but they are not - // convertible into each other, for example, and share no code (both - // instantiations are compiled completely independently). - - - // @sect4{LaplaceProblem::LaplaceProblem} - - // After this introduction, here is the constructor of the - // LaplaceProblem class. It specifies the desired polynomial - // degree of the finite elements and associates the DoFHandler to the - // triangulation just as in the previous example program, step-3: - template - LaplaceProblem::LaplaceProblem () : - fe (1), - dof_handler (triangulation) - {} - - - // @sect4{LaplaceProblem::make_grid_and_dofs} - - // Grid creation is something inherently dimension dependent. However, as - // long as the domains are sufficiently similar in 2D or 3D, the library can - // abstract for you. In our case, we would like to again solve on the square - // [-1,1]x[-1,1] in 2D, or on the cube [-1,1]x[-1,1]x[-1,1] in 3D; both can - // be termed hyper_cube, so we may use the same function in - // whatever dimension we are. Of course, the functions that create a - // hypercube in two and three dimensions are very much different, but that - // is something you need not care about. Let the library handle the - // difficult things. - // - // Likewise, associating a degree of freedom with each vertex is something - // which certainly looks different in 2D and 3D, but that does not need to - // bother you either. This function therefore looks exactly like in the - // previous example, although it performs actions that in their details are - // quite different if dim happens to be 3. The only significant - // difference from a user's perspective is the number of cells resulting, - // which is much higher in three than in two space dimensions! - template - void LaplaceProblem::make_grid_and_dofs () + template + void HeatEquation::run() { - GridGenerator::hyper_cube (triangulation, -30, 30); + setup_system(); - for (unsigned int f=0; f::faces_per_cell; ++f) - if (triangulation.begin()->face(f)->center()[2] > 15) - { - triangulation.begin()->face(f)->set_boundary_indicator (1); - for (unsigned int i=0; i::lines_per_face; ++i) - triangulation.begin()->face(f)->line(i)->set_boundary_indicator (1); - break; - } - triangulation.set_boundary (1, pds); + VectorTools::interpolate(dof_handler, ZeroFunction(), solution_u); + timestep_number = 0; + output_results(); - for (unsigned int v=0; v::vertices_per_cell; ++v) - if (triangulation.begin()->vertex(v)[2] > 0) - triangulation.begin()->vertex(v) - = pds.closest_point (Point<3>(triangulation.begin()->vertex(v)[0], - triangulation.begin()->vertex(v)[1], - 0)); + VectorTools::interpolate(dof_handler, ZeroFunction(), + old_solution_u); - for (unsigned int i=0; i<4; ++i) + Vector tmp(solution_u.size()); + Vector forcing_terms(solution_u.size()); + + for (timestep_number = 1, time = time_step; time <= 0.5; + time += time_step, ++timestep_number) { - for (typename Triangulation::active_cell_iterator - cell = triangulation.begin_active(); - cell != triangulation.end(); ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->boundary_indicator() == 1) - cell->set_refine_flag (); - - triangulation.execute_coarsening_and_refinement (); - - std::cout << "Refinement cycle " << i << std::endl - << " Number of active cells: " - << triangulation.n_active_cells() - << std::endl - << " Total number of cells: " - << triangulation.n_cells() + std::cout << "Time step " << timestep_number << " at t=" << time << std::endl; - } - - - dof_handler.distribute_dofs (fe); - - std::cout << " Number of degrees of freedom: " - << dof_handler.n_dofs() - << std::endl; - - 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); - sparsity_pattern.compress(); - - system_matrix.reinit (sparsity_pattern); - - solution.reinit (dof_handler.n_dofs()); - system_rhs.reinit (dof_handler.n_dofs()); - } - + mass_matrix.vmult(system_rhs, old_solution_u); - // @sect4{LaplaceProblem::assemble_system} - - // Unlike in the previous example, we would now like to use a non-constant - // right hand side function and non-zero boundary values. Both are tasks - // that are readily achieved with a only a few new lines of code in the - // assemblage of the matrix and right hand side. - // - // More interesting, though, is the way we assemble matrix and right hand - // side vector dimension independently: there is simply no difference to the - // two-dimensional case. Since the important objects used in this function - // (quadrature formula, FEValues) depend on the dimension by way of a - // template parameter as well, they can take care of setting up properly - // everything for the dimension for which this function is compiled. By - // declaring all classes which might depend on the dimension using a - // template parameter, the library can make nearly all work for you and you - // don't have to care about most things. - template - void LaplaceProblem::assemble_system () - { - MatrixTools::create_laplace_matrix (dof_handler, - QGauss(2), - system_matrix); - system_rhs = 0; - - std::map boundary_values; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - BoundaryValues(), - boundary_values); - MatrixTools::apply_boundary_values (boundary_values, - system_matrix, - solution, - system_rhs); - } + laplace_matrix.vmult(tmp, old_solution_u); + system_rhs.add(-(1 - theta) * time_step, tmp); //I omit here a time_step + RightHandSide rhs_function; + rhs_function.set_time(time); + VectorTools::create_right_hand_side(dof_handler, QGauss(2), + rhs_function, tmp); + forcing_terms = tmp; + forcing_terms *= theta; // I omit here a time_step - // @sect4{LaplaceProblem::solve} + rhs_function.set_time(time - time_step); + VectorTools::create_right_hand_side(dof_handler, QGauss(2), + rhs_function, tmp); - // Solving the linear system of equations is something that looks almost - // identical in most programs. In particular, it is dimension independent, - // so this function is copied verbatim from the previous example. - template - void LaplaceProblem::solve () - { - // NEW - SolverControl solver_control (dof_handler.n_dofs(), - 1e-12*system_rhs.l2_norm()); - SolverCG<> cg (solver_control); + forcing_terms.add((1 - theta) * time_step, tmp); - PreconditionSSOR<> preconditioner; - preconditioner.initialize(system_matrix, 1.2); - - cg.solve (system_matrix, solution, system_rhs, - preconditioner); - } - - - // @sect4{LaplaceProblem::output_results} - - // This function also does what the respective one did in step-3. No changes - // here for dimension independence either. - // - // The only difference to the previous example is that we want to write - // output in GMV format, rather than for gnuplot (GMV is another graphics - // program that, contrary to gnuplot, shows data in nice colors, allows - // rotation of geometries with the mouse, and generates reasonable - // representations of 3d data; for ways to obtain it see the ReadMe file of - // deal.II). To write data in this format, we simply replace the - // data_out.write_gnuplot call by - // data_out.write_gmv. - // - // Since the program will run both 2d and 3d versions of the laplace solver, - // we use the dimension in the filename to generate distinct filenames for - // each run (in a better program, one would check whether `dim' can have - // other values than 2 or 3, but we neglect this here for the sake of - // brevity). - template - void LaplaceProblem::output_results () const - { - DataOut data_out; - - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (solution, "solution"); - - data_out.build_patches (); - - std::ofstream output (dim == 2 ? - "solution-2d.gmv" : - "solution-3d.gmv"); - data_out.write_gmv (output); - } + system_rhs.add(time_step, forcing_terms); + { + BoundaryValuesU boundary_values_u_function; + boundary_values_u_function.set_time(time); + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, 0, + boundary_values_u_function, boundary_values); - // @sect4{LaplaceProblem::run} + matrix_u.copy_from(mass_matrix); + matrix_u.add(theta * time_step, laplace_matrix); // I omit here a time_step*theta + MatrixTools::apply_boundary_values(boundary_values, matrix_u, + solution_u, system_rhs); + } + solve_u(); - // This is the function which has the top-level control over - // everything. Apart from one line of additional output, it is the same as - // for the previous example. - template - void LaplaceProblem::run () - { - std::cout << "Solving problem in " << dim << " space dimensions." << std::endl; + output_results(); - make_grid_and_dofs(); - assemble_system (); - solve (); - output_results (); + old_solution_u = solution_u; + } } } - -// @sect3{The main function} - -// And this is the main function. It also looks mostly like in step-3, but if -// you look at the code below, note how we first create a variable of type -// LaplaceProblem@<2@> (forcing the compiler to compile the class -// template with dim replaced by 2) and run a 2d -// simulation, and then we do the whole thing over in 3d. -// -// In practice, this is probably not what you would do very frequently (you -// probably either want to solve a 2d problem, or one in 3d, but not both at -// the same time). However, it demonstrates the mechanism by which we can -// simply change which dimension we want in a single place, and thereby force -// the compiler to recompile the dimension independent class templates for the -// dimension we request. The emphasis here lies on the fact that we only need -// to change a single place. This makes it rather trivial to debug the program -// in 2d where computations are fast, and then switch a single place to a 3 to -// run the much more computing intensive program in 3d for `real' -// computations. -// -// Each of the two blocks is enclosed in braces to make sure that the -// laplace_problem_2d variable goes out of scope (and releases -// the memory it holds) before we move on to allocate memory for the 3d -// case. Without the additional braces, the laplace_problem_2d -// variable would only be destroyed at the end of the function, i.e. after -// running the 3d problem, and would needlessly hog memory while the 3d run -// could actually use it. -// -// Finally, the first line of the function is used to suppress some output. -// Remember that in the previous example, we had the output from the linear -// solvers about the starting residual and the number of the iteration where -// convergence was detected. This can be suppressed through the -// deallog.depth_console(0) call. -// -// The rationale here is the following: the deallog (i.e. deal-log, not -// de-allog) variable represents a stream to which some parts of the library -// write output. It redirects this output to the console and if required to a -// file. The output is nested in a way so that each function can use a prefix -// string (separated by colons) for each line of output; if it calls another -// function, that may also use its prefix which is then printed after the one -// of the calling function. Since output from functions which are nested deep -// below is usually not as important as top-level output, you can give the -// deallog variable a maximal depth of nested output for output to console and -// file. The depth zero which we gave here means that no output is written. By -// changing it you can get more information about the innards of the library. -int main () +int main() { try { using namespace dealii; using namespace Step26; - deallog.depth_console (0); + deallog.depth_console(0); + + HeatEquation<2> heat_equation_solver; + heat_equation_solver.run(); - LaplaceProblem<3> laplace_problem_3d; - laplace_problem_3d.run (); } catch (std::exception &exc) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl + std::cerr << "Exception on processing: " << std::endl << exc.what() + << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; @@ -677,8 +307,8 @@ int main () std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl + std::cerr << "Unknown exception!" << std::endl << "Aborting!" + << std::endl << "----------------------------------------------------" << std::endl; return 1; -- 2.39.5