<map id="StepsMap" name="StepsMap">
-<area shape="poly" id="node1" href="../deal.II/step_1.html" title="Creating a grid. Refining it. Writing it to a file" alt="" coords="548,60,548,28,525,5,493,5,471,28,471,60,493,83,525,83"/>
-<area shape="poly" id="node2" href="../deal.II/step_2.html" title="Assigning degrees of freedom to a grid." alt="" coords="548,185,548,153,525,131,493,131,471,153,471,185,493,208,525,208"/>
-<area shape="poly" id="node3" href="../deal.II/step_3.html" title="Solving Poisson's equation." alt="" coords="548,311,548,279,525,256,493,256,471,279,471,311,493,333,525,333"/>
-<area shape="poly" id="node4" href="../deal.II/step_4.html" title="Dimension independent programming. Boundary conditions." alt="" coords="548,436,548,404,525,381,493,381,471,404,471,436,493,459,525,459"/>
-<area shape="poly" id="node5" href="../deal.II/step_5.html" title="Reading a grid from disk. Computations on successively refined grids." alt="" coords="437,561,437,529,415,507,383,507,360,529,360,561,383,584,415,584"/>
-<area shape="rect" id="node10" href="../deal.II/step_10.html" title="Higher order mappings." alt="" coords="239,529,275,562"/>
-<area shape="rect" id="node15" href="../deal.II/step_15.html" title="1d problems. A nonlinear problem." alt="" coords="299,529,335,562"/>
-<area shape="rect" id="node20" href="../deal.II/step_20.html" title="Mixed finite elements for the mixed Laplacian. Block solvers." alt="" coords="462,533,497,558"/>
-<area shape="rect" id="node23" href="../deal.II/step_23.html" title="Time dependent problems. The wave equation." alt="" coords="146,533,181,558"/>
-<area shape="rect" id="node28" href="../deal.II/step_29.html" title="A complex-valued Helmholtz equation. Sparse direct solvers." alt="" coords="522,529,557,562"/>
-<area shape="rect" id="node33" href="../deal.II/step_34.html" title="Boundary element methods for potential flow." alt="" coords="582,533,617,558"/>
-<area shape="rect" id="node35" href="../deal.II/step_36.html" title="Finding eigenvalues of the Schrödinger equation." alt="" coords="642,529,677,562"/>
-<area shape="rect" id="node37" href="../deal.II/step_38.html" title="Solve the Laplace Beltrami operator on a Half Sphere." alt="" coords="625,654,660,687"/>
-<area shape="rect" id="node40" href="../deal.II/step_41.html" title="Solving the obstacle problem (a variational inequality)" alt="" coords="754,533,789,558"/>
-<area shape="poly" id="node6" href="../deal.II/step_6.html" title="Adaptive local refinement. Higher order elements" alt="" coords="408,687,408,655,385,632,353,632,331,655,331,687,353,709,385,709"/>
-<area shape="rect" id="node7" href="../deal.II/step_7.html" title="Helmholtz equation. Computing errors. Boundary integrals." alt="" coords="153,758,185,791"/>
-<area shape="rect" id="node8" href="../deal.II/step_8.html" title="Systems of PDE. Elasticity." alt="" coords="442,758,475,791"/>
-<area shape="rect" id="node9" href="../deal.II/step_9.html" title="Advection equation. Multithreading. Refinement criteria." alt="" coords="499,758,532,791"/>
-<area shape="rect" id="node13" href="../deal.II/step_13.html" title="Modularity. Software design." alt="" coords="41,758,76,791"/>
-<area shape="rect" id="node16" href="../deal.II/step_16.html" title="Multigrid on adaptive meshes." alt="" coords="101,841,136,874"/>
-<area shape="rect" id="node22" href="../deal.II/step_22.html" title="The Stokes equation on adaptive meshes." alt="" coords="609,762,644,787"/>
-<area shape="rect" id="node26" href="../deal.II/step_27.html" title="hp-adaptive finite element methods." alt="" coords="322,758,357,791"/>
-<area shape="rect" id="node27" href="../deal.II/step_28.html" title="Handling multiple meshes at the same time. Neutron transport." alt="" coords="210,758,245,791"/>
-<area shape="rect" id="node38" href="../deal.II/step_39.html" title="Interior Penalty for the Laplace equation. Adaptive refinement. Multigrid." alt="" coords="265,923,300,957"/>
-<area shape="rect" id="node39" href="../deal.II/step_40.html" title="Solving the Laplace equation on adaptive meshes on thousands of processors." alt="" coords="553,923,588,957"/>
-<area shape="rect" id="node43" href="../deal.II/step_45.html" title="Periodic boundary conditions" alt="" coords="382,758,417,791"/>
-<area shape="rect" id="node12" href="../deal.II/step_12.html" title="Discontinuous Galerkin for linear advection." alt="" coords="205,841,240,874"/>
-<area shape="rect" id="node17" href="../deal.II/step_17.html" title="Parallel computing using MPI. Using PETSc." alt="" coords="497,845,532,870"/>
-<area shape="rect" id="node42" href="../deal.II/step_44.html" title="Quasi-static finite-strain elasticity" alt="" coords="378,1010,413,1035"/>
-<area shape="rect" id="node44" href="../deal.II/step_46.html" title="Coupling different physical models (flow, elasticity) in different parts of the domain" alt="" coords="377,841,412,874"/>
-<area shape="rect" id="node11" href="../deal.II/step_11.html" title="Higher order mappings. Dealing with constraints." alt="" coords="239,654,275,687"/>
-<area shape="rect" id="node29" href="../deal.II/step_30.html" title="Anisotropic refinement for DG methods." alt="" coords="205,923,240,957"/>
-<area shape="rect" id="node32" href="../deal.II/step_33.html" title="Hyperbolic conservation laws: the Euler equations of gas dynamics." alt="" coords="145,927,180,953"/>
-<area shape="rect" id="node14" href="../deal.II/step_14.html" title="Duality based error estimates. Adaptivity." alt="" coords="41,841,76,874"/>
-<area shape="rect" id="node36" href="../deal.II/step_37.html" title="Matrix-free methods. Multigrid. Cell-based finite element operator." alt="" coords="85,923,120,957"/>
-<area shape="rect" id="node18" href="../deal.II/step_18.html" title="Quasistatic elasticity. More parallel computing." alt="" coords="437,927,472,953"/>
-<area shape="rect" id="node31" href="../deal.II/step_32.html" title="A parallel Boussinesq flow solver for thermal convection in the earth mantle." alt="" coords="553,1010,588,1035"/>
-<area shape="rect" id="node19" href="../deal.II/step_19.html" title="Handling input parameter files. Converting output formats." alt="" coords="437,841,472,874"/>
-<area shape="rect" id="node21" href="../deal.II/step_21.html" title="Two-phase flow in porous media." alt="" coords="565,658,600,683"/>
-<area shape="rect" id="node41" href="../deal.II/step_43.html" title="Efficient ways to solve two-phase flow problems on adaptive meshes in 2d and 3d." alt="" coords="697,927,732,953"/>
-<area shape="rect" id="node30" href="../deal.II/step_31.html" title="Boussinesq flow for thermal convection." alt="" coords="669,845,704,870"/>
-<area shape="rect" id="node34" href="../deal.II/step_35.html" title="A projection solver for the Navier-Stokes equations." alt="" coords="609,845,644,870"/>
-<area shape="rect" id="node24" href="../deal.II/step_24.html" title="The wave equation with absorbing boundary conditions. Extracting point values." alt="" coords="53,658,88,683"/>
-<area shape="rect" id="node25" href="../deal.II/step_25.html" title="The nonlinear sine-Gordon soliton equation" alt="" coords="5,927,40,953"/>
-<area shape="rect" id="node45" href="../deal.II/step_48.html" title="Parallelization via MPI. The wave equation, in linear and nonlinear variants. Mass lumping. Cell-based finite element operator." alt="" coords="85,1006,120,1039"/>
+<area shape="poly" id="node1" href="../deal.II/step_1.html" title="Creating a grid. Refining it. Writing it to a file" alt="" coords="572,60 572,28 549,5 517,5 495,28 495,60 517,83 549,83"/>
+<area shape="poly" id="node2" href="../deal.II/step_2.html" title="Assigning degrees of freedom to a grid." alt="" coords="572,185 572,153 549,131 517,131 495,153 495,185 517,208 549,208"/>
+<area shape="poly" id="node3" href="../deal.II/step_3.html" title="Solving Poisson's equation." alt="" coords="572,311 572,279 549,256 517,256 495,279 495,311 517,333 549,333"/>
+<area shape="poly" id="node4" href="../deal.II/step_4.html" title="Dimension independent programming. Boundary conditions." alt="" coords="572,436 572,404 549,381 517,381 495,404 495,436 517,459 549,459"/>
+<area shape="poly" id="node5" href="../deal.II/step_5.html" title="Reading a grid from disk. Computations on successively refined grids." alt="" coords="404,561 404,529 381,507 349,507 327,529 327,561 349,584 381,584"/>
+<area shape="rect" id="node10" href="../deal.II/step_10.html" title="Higher order mappings." alt="" coords="487,529,521,562"/>
+<area shape="rect" id="node20" href="../deal.II/step_20.html" title="Mixed finite elements for the mixed Laplacian. Block solvers." alt="" coords="164,531,199,560"/>
+<area shape="rect" id="node23" href="../deal.II/step_23.html" title="Time dependent problems. The wave equation." alt="" coords="428,531,463,560"/>
+<area shape="rect" id="node28" href="../deal.II/step_29.html" title="A complex-valued Helmholtz equation. Sparse direct solvers." alt="" coords="545,529,580,562"/>
+<area shape="rect" id="node33" href="../deal.II/step_34.html" title="Boundary element methods for potential flow." alt="" coords="655,531,689,560"/>
+<area shape="rect" id="node35" href="../deal.II/step_36.html" title="Finding eigenvalues of the Schrödinger equation." alt="" coords="713,529,748,562"/>
+<area shape="rect" id="node37" href="../deal.II/step_38.html" title="Solve the Laplace Beltrami operator on a Half Sphere." alt="" coords="655,654,689,687"/>
+<area shape="poly" id="node6" href="../deal.II/step_6.html" title="Adaptive local refinement. Higher order elements" alt="" coords="404,687 404,655 381,632 349,632 327,655 327,687 349,709 381,709"/>
+<area shape="rect" id="node7" href="../deal.II/step_7.html" title="Helmholtz equation. Computing errors. Boundary integrals." alt="" coords="549,758,581,791"/>
+<area shape="rect" id="node8" href="../deal.II/step_8.html" title="Systems of PDE. Elasticity." alt="" coords="151,758,183,791"/>
+<area shape="rect" id="node9" href="../deal.II/step_9.html" title="Advection equation. Multithreading. Refinement criteria." alt="" coords="207,758,239,791"/>
+<area shape="rect" id="node13" href="../deal.II/step_13.html" title="Modularity. Software design." alt="" coords="431,758,465,791"/>
+<area shape="rect" id="node15" href="../deal.II/step_15.html" title="A nonlinear elliptic problem. Newton's method." alt="" coords="489,758,524,791"/>
+<area shape="rect" id="node16" href="../deal.II/step_16.html" title="Multigrid on adaptive meshes." alt="" coords="657,841,692,874"/>
+<area shape="rect" id="node22" href="../deal.II/step_22.html" title="The Stokes equation on adaptive meshes." alt="" coords="33,760,68,789"/>
+<area shape="rect" id="node26" href="../deal.II/step_27.html" title="hp-adaptive finite element methods." alt="" coords="92,758,127,791"/>
+<area shape="rect" id="node27" href="../deal.II/step_28.html" title="Handling multiple meshes at the same time. Neutron transport." alt="" coords="263,758,297,791"/>
+<area shape="rect" id="node38" href="../deal.II/step_39.html" title="Interior Penalty for the Laplace equation. Adaptive refinement. Multigrid." alt="" coords="597,923,632,957"/>
+<area shape="rect" id="node39" href="../deal.II/step_40.html" title="Solving the Laplace equation on adaptive meshes on thousands of processors." alt="" coords="376,923,411,957"/>
+<area shape="rect" id="node43" href="../deal.II/step_45.html" title="Periodic boundary conditions" alt="" coords="321,758,356,791"/>
+<area shape="rect" id="node12" href="../deal.II/step_12.html" title="Discontinuous Galerkin for linear advection." alt="" coords="548,841,583,874"/>
+<area shape="rect" id="node17" href="../deal.II/step_17.html" title="Parallel computing using MPI. Using PETSc." alt="" coords="267,843,301,872"/>
+<area shape="rect" id="node42" href="../deal.II/step_44.html" title="Quasi-static finite-strain elasticity" alt="" coords="296,1005,331,1035"/>
+<area shape="rect" id="node44" href="../deal.II/step_46.html" title="Coupling different physical models (flow, elasticity) in different parts of the domain" alt="" coords="149,841,184,874"/>
+<area shape="rect" id="node11" href="../deal.II/step_11.html" title="Higher order mappings. Dealing with constraints." alt="" coords="487,654,521,687"/>
+<area shape="rect" id="node29" href="../deal.II/step_30.html" title="Anisotropic refinement for DG methods." alt="" coords="480,923,515,957"/>
+<area shape="rect" id="node32" href="../deal.II/step_33.html" title="Hyperbolic conservation laws: the Euler equations of gas dynamics." alt="" coords="539,925,573,955"/>
+<area shape="rect" id="node14" href="../deal.II/step_14.html" title="Duality based error estimates. Adaptivity." alt="" coords="431,841,465,874"/>
+<area shape="rect" id="node40" href="../deal.II/step_41.html" title="Solving the obstacle problem (a variational inequality)" alt="" coords="489,843,524,872"/>
+<area shape="rect" id="node36" href="../deal.II/step_37.html" title="Matrix-free methods. Multigrid. Cell-based finite element operator." alt="" coords="656,923,691,957"/>
+<area shape="rect" id="node18" href="../deal.II/step_18.html" title="Quasistatic elasticity. More parallel computing." alt="" coords="267,925,301,955"/>
+<area shape="rect" id="node31" href="../deal.II/step_32.html" title="A parallel Boussinesq flow solver for thermal convection in the earth mantle." alt="" coords="212,1005,247,1035"/>
+<area shape="rect" id="node19" href="../deal.II/step_19.html" title="Handling input parameter files. Converting output formats." alt="" coords="208,841,243,874"/>
+<area shape="rect" id="node21" href="../deal.II/step_21.html" title="Two-phase flow in porous media." alt="" coords="33,656,68,685"/>
+<area shape="rect" id="node41" href="../deal.II/step_43.html" title="Efficient ways to solve two-phase flow problems on adaptive meshes in 2d and 3d." alt="" coords="5,925,40,955"/>
+<area shape="rect" id="node30" href="../deal.II/step_31.html" title="Boussinesq flow for thermal convection." alt="" coords="32,843,67,872"/>
+<area shape="rect" id="node34" href="../deal.II/step_35.html" title="A projection solver for the Navier-Stokes equations." alt="" coords="91,843,125,872"/>
+<area shape="rect" id="node24" href="../deal.II/step_24.html" title="The wave equation with absorbing boundary conditions. Extracting point values." alt="" coords="428,656,463,685"/>
+<area shape="rect" id="node25" href="../deal.II/step_25.html" title="The nonlinear sine-Gordon soliton equation" alt="" coords="715,925,749,955"/>
+<area shape="rect" id="node45" href="../deal.II/step_48.html" title="Parallelization via MPI. The wave equation, in linear and nonlinear variants. Mass lumping. Cell-based finite element operator." alt="" coords="656,1005,691,1035"/>
</map>
#include <fstream>
#include <iostream>
- // We will use adaptive mesh refinement between Newton
- // interations. To do so, we need to be able to work
- // with a solution on the new mesh, although it was
- // computed on the old one. The SolutionTransfer
- // class transfers the solution to the new mesh.
+ // We will use adaptive mesh refinement
+ // between Newton interations. To do so, we
+ // need to be able to work with a solution on
+ // the new mesh, although it was computed on
+ // the old one. The SolutionTransfer class
+ // transfers the solution from the old to the
+ // new mesh:
#include <deal.II/numerics/solution_transfer.h>
- // As in previous programs:
+ // We then open a namepsace for this program
+ // and import everything from the dealii
+ // namespace into it, as in previous
+ // programs:
namespace Step15
{
using namespace dealii;
// @sect3{The <code>MinimalSurfaceProblem</code> class template}
- // The class template is basically
- // the same as in step 6. Four
- // additions are made: There are
- // two solution vectors, one for
- // the Newton update, and one for
- // the solution of the original
- // pde. Also we need a double for
- // the residual of the Newton
- // method, an integer, which counts
- // the mesh refinements and a bool
- // for the boundary condition in
- // the first Newton step.
+ // The class template is basically the same
+ // as in step-6. Four additions are made:
+ // - There are two solution vectors, one for
+ // the Newton update $\delta u^n$, and one
+ // for the current iterate $u^n$.
+ // - The <code>setup_system</code> function
+ // takes an argument that denotes whether
+ // this is the first time it is called or
+ // not. The difference is that the first
+ // time around we need to distributed
+ // degrees of freedom and set the
+ // solution vector for $u^n$ to the
+ // correct size. The following times, the
+ // function is called after we have
+ // already done these steps as part of
+ // refining the mesh in
+ // <code>refine_mesh</code>.
+ // - We then also need new functions:
+ // <code>set_boundary_values()</code>
+ // takes care of setting the boundary
+ // values on the solution vector
+ // correctly, as discussed at the end of
+ // the
+ // introduction. <code>compute_residual()</code>
+ // is a function that computes the norm
+ // of the nonlinear (discrete)
+ // residual. We use this function to
+ // monitor convergence of the Newton
+ // iteration. The function takes a step
+ // length $\alpha^n$ as argument to
+ // compute the residual of $u^n +
+ // \alpha^n \; \delta u^n$. This is
+ // something one typically needs for step
+ // length control, although we will not
+ // use this feature here. Finally,
+ // <code>determine_step_length()</code>
+ // computes the step length $\alpha^n$ in
+ // each Newton iteration. As discussed in
+ // the introduction, we here use a fixed
+ // step length and leave implementing a
+ // better strategy as an exercise.
template <int dim>
class MinimalSurfaceProblem
void run ();
private:
- void setup_system ();
+ void setup_system (const bool initial_step);
void assemble_system ();
void solve ();
- void refine_grid ();
+ void refine_mesh ();
void set_boundary_values ();
double compute_residual (const double alpha) const;
double determine_step_length () const;
Vector<double> present_solution;
Vector<double> newton_update;
Vector<double> system_rhs;
-
-
-
- unsigned int refinement;
-
- // As described in the
- // Introduction, the first
- // Newton iteration is special,
- // because of the boundary
- // condition. To implement
- // these correctly, there is a
- // bool, which is true in the
- // first step and false ever
- // after.
- bool first_step;
};
// @sect3{Boundary condition}
// The boundary condition is
- // implemented just like in step 4.
- // It was chosen as $g(x,y)=sin(2
- // \pi (x+y))$ in this example.
+ // implemented just like in step-4.
+ // It is chosen as $g(x,y)=\sin(2
+ // \pi (x+y))$:
template <int dim>
class BoundaryValues : public Function<dim>
const unsigned int component = 0) const;
};
+
template <int dim>
double BoundaryValues<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
- //
template <int dim>
MinimalSurfaceProblem<dim>::~MinimalSurfaceProblem ()
{
// @sect4{MinimalSurfaceProblem::setup_system}
- // As always in the setup-system
- // function, we setup the variables
- // of the finite element
- // method. There are same
- // differences to step 6, because
- // we don't have to solve one pde
- // over all, but one in every
- // Newton step. Also the starting
- // function has to be setup in the
- // first step.
+ // As always in the setup-system function,
+ // we setup the variables of the finite
+ // element method. There are same
+ // differences to step-6, because there we
+ // start solving the PDE from scratch in
+ // every refinement cycle whereas here we
+ // need to take the solution from the
+ // previous mesh onto the current
+ // mesh. Consequently, we can't just reset
+ // solution vectors. The argument passed to
+ // this function thus indicates whether we
+ // can distributed degrees of freedom (plus
+ // compute constraints) and set the
+ // solution vector to zero or whether this
+ // has happened elsewhere already
+ // (specifically, in
+ // <code>refine_mesh()</code>).
template <int dim>
- void MinimalSurfaceProblem<dim>::setup_system ()
+ void MinimalSurfaceProblem<dim>::setup_system (const bool initial_step)
{
-
- // This function will be called,
- // every time we refine the mesh
- // to resize the system matrix,
- // Newton update - and right hand
- // side vector and to set the
- // right values of hanging nodes
- // to get a continuous solution.
- // But only the first time, the
- // starting solution has to be
- // initialized. Also the vector
- // of the solution will be
- // resized in the
- // <code>refine_grid</code>
- // function, while the vector is
- // transferred to the new mesh.
-
- if (first_step)
+ if (initial_step)
{
dof_handler.distribute_dofs (fe);
present_solution.reinit (dof_handler.n_dofs());
- // The constraint matrix,
- // holding a list of the
- // hanging nodes, will be
- // setup in the
- // <code>refine_grid</code>
- // function after refining
- // the mesh.
-
hanging_node_constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler,
hanging_node_constraints);
// The remaining parts of the
// function are the same as in
- // step 6.
+ // step-6.
newton_update.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
// @sect4{MinimalSurfaceProblem::assemble_system}
- // This function does the same as
- // in the previous tutorials. The
- // only additional step is the
- // correct implementation of the
- // boundary condition and the usage
- // of the gradients of the old
- // solution.
-
+ // This function does the same as in the
+ // previous tutorials except that now, of
+ // course, the matrix and right hand side
+ // functions depend on the previous
+ // iteration's solution. As discussed in
+ // the introduction, we need to use zero
+ // boundary values for the Newton updates;
+ // we compute them at the end of this
+ // function.
+ //
+ // The top of the function contains the
+ // usual boilerplate code, setting up the
+ // objects that allow us to evaluate shape
+ // functions at quadrature points and
+ // temporary storage locations for the
+ // local matrices and vectors, as well as
+ // for the gradients of the previous
+ // solution at the quadrature points. We
+ // then start the loop over all cells:
template <int dim>
void MinimalSurfaceProblem<dim>::assemble_system ()
{
system_rhs = 0;
FEValues<dim> fe_values (fe, quadrature_formula,
- update_gradients |
- update_quadrature_points | update_JxW_values);
+ update_gradients |
+ update_quadrature_points |
+ update_JxW_values);
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- const unsigned int n_q_points = quadrature_formula.size();
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- Vector<double> cell_rhs (dofs_per_cell);
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
- std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ std::vector<Tensor<1, dim> > old_solution_gradients(n_q_points);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
fe_values.reinit (cell);
-
+ // For the assembly of the linear
+ // system, we have to obtain the
+ // values of the previous solution's
+ // gradients at the quadrature
+ // points. There is a standard way of
+ // doing this: the
+ // FEValues::get_function function
+ // takes a vector that represents a
+ // finite element field defined on a
+ // DoFHandler, and evaluates the
+ // gradients of this field at the
+ // quadrature points of the cell with
+ // which the FEValues object has last
+ // been reinitialized. The values of
+ // the gradients at all quadrature
+ // points are then written into the
+ // second argument:
+ fe_values.get_function_gradients(present_solution,
+ old_solution_gradients);
+
+ // With this, we can then do the
+ // integration loop over all
+ // quadrature points and shape
+ // functions. Having just computed
+ // the gradients of the old solution
+ // in the quadrature points, we are
+ // able to compute the coefficients
+ // $a_{n}$ in these points. The
+ // assembly of the system itself then
+ // looks similar to what we always do
+ // with the exception of the
+ // nonlinear terms, as does copying
+ // the results from the local objects
+ // into the global ones:
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
+ const double coeff
+ = 1.0 / std::sqrt(1 +
+ old_solution_gradients[q_point] *
+ old_solution_gradients[q_point]);
- // To setup up the linear
- // system, the gradient of
- // the old solution in the
- // quadrature points is
- // needed. For this purpose
- // there is is a function,
- // which will write these
- // gradients in a vector,
- // where every component of
- // the vector is a vector
- // itself:
-
- std::vector<Tensor<1, dim> > gradients(n_q_points);
- fe_values.get_function_gradients(present_solution, gradients);
-
- // Having the gradients of
- // the old solution in the
- // quadrature points, we
- // are able to compute the
- // coefficients $a_{n}$ in
- // these points.
-
- const double coeff = 1/sqrt(1 + gradients[q_point] * gradients[q_point]);
-
- // The assembly of the
- // system then is the same
- // as always, except of the
- // damping parameter of the
- // Newton method, which we
- // set on 0.1 in this case.
-
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
{
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
{
cell_matrix(i, j) += (fe_values.shape_grad(i, q_point)
* coeff
* (fe_values.shape_grad(j, q_point)
- - coeff * coeff
+ -
+ coeff * coeff
* (fe_values.shape_grad(j, q_point)
- * gradients[q_point])
- * gradients[q_point])
+ *
+ old_solution_gradients[q_point])
+ * old_solution_gradients[q_point]
+ )
* fe_values.JxW(q_point));
}
- cell_rhs(i) -= (fe_values.shape_grad(i, q_point) * coeff
- * gradients[q_point] * fe_values.JxW(q_point));
+ cell_rhs(i) -= (fe_values.shape_grad(i, q_point)
+ * coeff
+ * old_solution_gradients[q_point]
+ * fe_values.JxW(q_point));
}
}
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
+
+ // Finally, we remove hanging nodes from
+ // the system and apply zero boundary
+ // values to the linear system that
+ // defines the Newton updates $\delta
+ // u^n$:
hanging_node_constraints.condense (system_matrix);
hanging_node_constraints.condense (system_rhs);
- std::map<unsigned int,double> boundary_values;
-
+ std::map<unsigned int,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
ZeroFunction<dim>(),
boundary_values);
-
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
newton_update,
}
- template <int dim>
- double MinimalSurfaceProblem<dim>::compute_residual (const double alpha) const
- {
- const QGauss<dim> quadrature_formula(3);
-
- Vector<double> residual (dof_handler.n_dofs());
-
- Vector<double> linearization_point (dof_handler.n_dofs());
- linearization_point = present_solution;
- linearization_point.add (alpha, newton_update);
-
- FEValues<dim> fe_values (fe, quadrature_formula,
- update_gradients |
- update_quadrature_points | update_JxW_values);
-
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- const unsigned int n_q_points = quadrature_formula.size();
-
- Vector<double> cell_rhs (dofs_per_cell);
-
- std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
- typename DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell!=endc; ++cell)
- {
- cell_rhs = 0;
-
- fe_values.reinit (cell);
-
-
- for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
- {
-
- // To setup up the linear
- // system, the gradient of
- // the old solution in the
- // quadrature points is
- // needed. For this purpose
- // there is is a function,
- // which will write these
- // gradients in a vector,
- // where every component of
- // the vector is a vector
- // itself:
-
- std::vector<Tensor<1, dim> > gradients(n_q_points);
- fe_values.get_function_gradients(linearization_point, gradients);
-
- // Having the gradients of
- // the old solution in the
- // quadrature points, we
- // are able to compute the
- // coefficients $a_{n}$ in
- // these points.
-
- const double coeff = 1/sqrt(1 + gradients[q_point] * gradients[q_point]);
-
- // The assembly of the
- // system then is the same
- // as always, except of the
- // damping parameter of the
- // Newton method, which we
- // set on 0.1 in this case.
-
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- cell_rhs(i) -= (fe_values.shape_grad(i, q_point) * coeff
- * gradients[q_point] * fe_values.JxW(q_point));
- }
-
- cell->get_dof_indices (local_dof_indices);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- residual(local_dof_indices[i]) += cell_rhs(i);
- }
- hanging_node_constraints.condense (residual);
-
- std::map<unsigned int,double> boundary_values;
- VectorTools::interpolate_boundary_values (dof_handler,
- 0,
- ZeroFunction<dim>(),
- boundary_values);
- for (std::map<unsigned int,double>::const_iterator p = boundary_values.begin();
- p != boundary_values.end(); ++p)
- residual(p->first) = 0;
-
- return residual.l2_norm();
- }
// @sect4{MinimalSurfaceProblem::solve}
- // The solve function is the same
- // as always, we just have to
- // implement the minimal residual
- // method as a solver and apply the
- // Newton update to the solution.
-
+ // The solve function is the same as
+ // always. At the end of the solution
+ // process we update the current solution
+ // by setting $u^{n+1}=u^n+\alpha^n\;\delta
+ // u^n$.
template <int dim>
void MinimalSurfaceProblem<dim>::solve ()
{
hanging_node_constraints.distribute (newton_update);
- // In this step, the old solution
- // is updated to the new one:
const double alpha = determine_step_length();
- std::cout << " step length alpha=" << alpha << std::endl;
present_solution.add (alpha, newton_update);
}
- template <int dim>
- double MinimalSurfaceProblem<dim>::determine_step_length() const
- {
- return 0.1;
- }
- // @sect4{MinimalSurfaceProblem::refine_grid}
-
- // The first part of this function
- // is the same as in step 6. But
- // after refining the mesh we have
- // to transfer the old solution to
- // the new one, which is done with
- // the help of the SolutionTransfer
- // class.
-
+ // @sect4{MinimalSurfaceProblem::refine_mesh}
+ // The first part of this function is the
+ // same as in step-6... However, after
+ // refining the mesh we have to transfer
+ // the old solution to the new one which we
+ // do with the help of the SolutionTransfer
+ // class. The process is slightly
+ // convoluted, so let us describe it in
+ // detail:
template <int dim>
- void MinimalSurfaceProblem<dim>::refine_grid ()
+ void MinimalSurfaceProblem<dim>::refine_mesh ()
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
estimated_error_per_cell,
0.3, 0.03);
- // Then we need an additional
- // step: if, for example, you
- // flag a cell that is once more
- // refined than its neighbor, and
- // that neighbor is not flagged
- // for refinement, we would end
- // up with a jump of two
- // refinement levels across a
- // cell interface. To avoid
- // these situations, the library
- // will silently also have to
- // refine the neighbor cell
- // once. It does so by calling
- // the
+ // Then we need an additional step: if,
+ // for example, you flag a cell that is
+ // once more refined than its neighbor,
+ // and that neighbor is not flagged for
+ // refinement, we would end up with a
+ // jump of two refinement levels across a
+ // cell interface. To avoid these
+ // situations, the library will silently
+ // also have to refine the neighbor cell
+ // once. It does so by calling the
// Triangulation::prepare_coarsening_and_refinement
- // function before actually doing
- // the refinement and coarsening.
- // This function flags a set of
- // additional cells for
- // refinement or coarsening, to
+ // function before actually doing the
+ // refinement and coarsening. This
+ // function flags a set of additional
+ // cells for refinement or coarsening, to
// enforce rules like the
- // one-hanging-node rule. The
- // cells that are flagged for
- // refinement and coarsening
- // after calling this function
+ // one-hanging-node rule. The cells that
+ // are flagged for refinement and
+ // coarsening after calling this function
// are exactly the ones that will
// actually be refined or
- // coarsened. Since the
- // SolutionTransfer class needs
- // this information in order to
- // store the data from the old
- // mesh and transfer to the new
- // one.
-
+ // coarsened. Usually, you don't have to
+ // do this by hand
+ // (Triangulation::execute_coarsening_and_refinement
+ // does this for you). However, we need
+ // to initialize the SolutionTransfer
+ // class and it needs to know the final
+ // set of cells that will be coarsened or
+ // refined in order to store the data
+ // from the old mesh and transfer to the
+ // new one. Thus, we call the function by
+ // hand:
triangulation.prepare_coarsening_and_refinement ();
// With this out of the way, we
- // initialize a SolutionTransfer
- // object with the present
- // DoFHandler and attach the
- // solution vector to it:
-
+ // initialize a SolutionTransfer object
+ // with the present DoFHandler and attach
+ // the solution vector to it, followed by
+ // doing the actual refinement and
+ // distribution of degrees of freedom on
+ // the new mesh
SolutionTransfer<dim> solution_transfer(dof_handler);
solution_transfer.prepare_for_coarsening_and_refinement(present_solution);
- // Then we do the actual
- // refinement, and distribute
- // degrees of freedom on the new
- // mesh:
-
triangulation.execute_coarsening_and_refinement();
- dof_handler.distribute_dofs(fe);
- // Finally, we retrieve the old
- // solution interpolated to the
- // new mesh. Since the
- // SolutionTransfer function does
- // not actually store the values
- // of the old solution, but
- // rather indices, we need to
- // preserve the old solution
- // vector until we have gotten
- // the new interpolated
- // values. Thus, we have the new
- // values written into a
- // temporary vector, and only
- // afterwards write them into the
- // solution vector object:
+ dof_handler.distribute_dofs(fe);
+ // Finally, we retrieve the old solution
+ // interpolated to the new mesh. Since
+ // the SolutionTransfer function does not
+ // actually store the values of the old
+ // solution, but rather indices, we need
+ // to preserve the old solution vector
+ // until we have gotten the new
+ // interpolated values. Thus, we have the
+ // new values written into a temporary
+ // vector, and only afterwards write them
+ // into the solution vector object. Once
+ // we have this solution we have to make
+ // sure that the $u^n$ we now have
+ // actually has the correct boundary
+ // values. As explained at the end of the
+ // introduction, this is not
+ // automatically the case even if the
+ // solution before refinement had the
+ // correct boundary values, and so we
+ // have to explicitly make sure that it
+ // now has:
Vector<double> tmp(dof_handler.n_dofs());
- solution_transfer.interpolate(present_solution,tmp);
- present_solution=tmp;
+ solution_transfer.interpolate(present_solution, tmp);
+ present_solution = tmp;
set_boundary_values ();
- // On the new mesh, there are
- // different hanging nodes, which
- // shall be enlisted in a matrix
- // like before. To ensure there
- // are no hanging nodes of the
- // old mesh in the matrix, it's
- // first cleared:
+ // On the new mesh, there are different
+ // hanging nodes, which we have to
+ // compute again. To ensure there are no
+ // hanging nodes of the old mesh in the
+ // object, it's first cleared. To be on
+ // the safe side, we then also make sure
+ // that the current solution's vector
+ // entries satisfy the hanging node
+ // constraints:
+
hanging_node_constraints.clear();
- // After doing so, the hanging
- // nodes of the new mesh can be
- // enlisted in the matrix, like
- // before. Calling the
- // <code>setup_system</code>
- // function in the
- // <code>run</code> function
- // again after this, the hanging
- // nodes don't have to be
- // enlisted there once more.
-
- DoFTools::make_hanging_node_constraints(dof_handler, hanging_node_constraints);
+ DoFTools::make_hanging_node_constraints(dof_handler,
+ hanging_node_constraints);
hanging_node_constraints.close();
- hanging_node_constraints.distribute(present_solution);
+
+ hanging_node_constraints.distribute (present_solution);
+
+ // We end the function by updating all
+ // the remaining data structures,
+ // indicating to
+ // <code>setup_dofs()</code> that this is
+ // not the first go-around and that it
+ // needs to preserve the content of the
+ // solution vector:
+ setup_system (false);
}
+
+ // @sect4{MinimalSurfaceProblem::set_boundary_values}
+
+ // The next function ensures that the
+ // solution vector's entries respect the
+ // boundary values for our problem. Having
+ // refined the mesh (or just started
+ // computations), there might be new nodal
+ // points on the boundary. These have
+ // values that are simply interpolated from
+ // the previous mesh (or are just zero),
+ // instead of the correct boundary
+ // values. This is fixed up by setting all
+ // boundary nodes explicit to the right
+ // value:
template <int dim>
void MinimalSurfaceProblem<dim>::set_boundary_values ()
{
- // Having refined the mesh, there
- // might be new nodal points on
- // the boundary. These have just
- // interpolated values, but not
- // the right boundary
- // values. This is fixed up, by
- // setting all boundary nodals
- // explicit to the right value:
-
- std::map<unsigned int, double> boundary_values2;
- VectorTools::interpolate_boundary_values(dof_handler, 0,
- BoundaryValues<dim>(), boundary_values2);
- for (std::map<unsigned int, double>::const_iterator p =
- boundary_values2.begin(); p != boundary_values2.end(); ++p)
+ std::map<unsigned int, double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ BoundaryValues<dim>(),
+ boundary_values);
+ for (std::map<unsigned int, double>::const_iterator
+ p = boundary_values.begin();
+ p != boundary_values.end(); ++p)
present_solution(p->first) = p->second;
}
- // @sect4{MinimalSurfaceProblem::run}
- // In the run function, the first
- // grid is build. Also in this
- // function, the Newton iteration
- // is implemented.
+ // @sect4{MinimalSurfaceProblem::compute_residual}
+
+ // In order to monitor convergence, we need
+ // a way to compute the norm of the
+ // (discrete) residual, i.e., the norm of
+ // the vector
+ // $\left<F(u^n),\varphi_i\right>$ with
+ // $F(u)=-\nabla \cdot \left(
+ // \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla
+ // u \right)$ as discussed in the
+ // introduction. It turns out that
+ // (although we don't use this feature in
+ // the current version of the program) one
+ // needs to compute the residual
+ // $\left<F(u^n+\alpha^n\;\delta u^n),\varphi_i\right>$
+ // when determining optimal step lengths,
+ // and so this is what we implement here:
+ // the function takes the step length
+ // $\alpha^n$ as an argument. The original
+ // functionality is of course obtained by
+ // passing a zero as argument.
+ //
+ // In the function below, we first set up a
+ // vector for the residual, and then a
+ // vector for the evaluation point
+ // $u^n+\alpha^n\;\delta u^n$. This is
+ // followed by the same boilerplate code we
+ // use for all integration operations:
template <int dim>
- void MinimalSurfaceProblem<dim>::run ()
+ double MinimalSurfaceProblem<dim>::compute_residual (const double alpha) const
{
+ Vector<double> residual (dof_handler.n_dofs());
- // The integer refinement counts
- // the mesh
- // refinements. Obviously
- // starting the program, it
- // should be zero.
- refinement=0;
- first_step=true;
-
- // As described in the
- // introduction, the domain is a
- // unitball around the
- // origin. The Mesh is globally
- // refined two times, not to
- // start on the coarse mesh,
- // which consists only of five
- // cells.
+ Vector<double> evaluation_point (dof_handler.n_dofs());
+ evaluation_point = present_solution;
+ evaluation_point.add (alpha, newton_update);
- GridGenerator::hyper_ball (triangulation);
- static const HyperBallBoundary<dim> boundary;
- triangulation.set_boundary (0, boundary);
- triangulation.refine_global(2);
+ const QGauss<dim> quadrature_formula(3);
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_gradients |
+ update_quadrature_points |
+ update_JxW_values);
- // The Newton iteration starts
- // here. During the first step,
- // there is no residual computed,
- // so the bool is needed here to
- // enter the iteration
- // scheme. Later the Newton
- // method will continue until the
- // residual is less than
- // $10^{-3}$.
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
- double previous_res = 0;
- while (first_step || (previous_res>1e-3))
+ Vector<double> cell_rhs (dofs_per_cell);
+ std::vector<Tensor<1, dim> > gradients(n_q_points);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
{
+ cell_rhs = 0;
+ fe_values.reinit (cell);
+
+ // The actual computation is much as
+ // in
+ // <code>assemble_system()</code>. We
+ // first evaluate the gradients of
+ // $u^n+\alpha^n\,\delta u^n$ at the
+ // quadrature points, then compute
+ // the coefficient $a_n$, and then
+ // plug it all into the formula for
+ // the residual:
+ fe_values.get_function_gradients (evaluation_point,
+ gradients);
- // In the first step, we
- // compute the solution on
- // the two times globally
- // refined mesh. After that
- // the mesh will be refined
- // adaptively, in order to
- // not get too many
- // cells. The refinement is
- // the first thing done every
- // time we restart the
- // process in the while-loop.
- if (!first_step)
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
{
- refine_grid();
+ const double coeff = 1/std::sqrt(1 +
+ gradients[q_point] *
+ gradients[q_point]);
- std::cout<<"********mesh-refinement:"<<refinement+1<<" ********"<<std::endl;
- refinement++;
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ cell_rhs(i) -= (fe_values.shape_grad(i, q_point)
+ * coeff
+ * gradients[q_point]
+ * fe_values.JxW(q_point));
}
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ residual(local_dof_indices[i]) += cell_rhs(i);
+ }
+
+ // At the end of this function we also
+ // have to deal with the hanging node
+ // constraints and with the issue of
+ // boundary values. With regard to the
+ // latter, we have to set to zero the
+ // elements of the residual vector for
+ // all entries that correspond to degrees
+ // of freedom that sit at the
+ // boundary. The reason is that because
+ // the value of the solution there is
+ // fixed, they are of course no "real"
+ // degrees of freedom and so, strictly
+ // speaking, we shouldn't have assembled
+ // entries in the residual vector for
+ // them. However, as we always do, we
+ // want to do exactly the same thing on
+ // every cell and so we didn't not want
+ // to deal with the question of whether a
+ // particular degree of freedom sits at
+ // the boundary in the integration
+ // above. Rather, we will simply set to
+ // zero these entries after the fact. To
+ // this end, we first need to determine
+ // which degrees of freedom do in fact
+ // belong to the boundary and then loop
+ // over all of those and set the residual
+ // entry to zero. This happens in the
+ // following lines:
+ hanging_node_constraints.condense (residual);
+
+ std::vector<bool> boundary_dofs (dof_handler.n_dofs());
+ DoFTools::extract_boundary_dofs (dof_handler,
+ std::vector<bool>(1,true),
+ boundary_dofs);
+ for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+ if (boundary_dofs[i] == true)
+ residual(i) = 0;
+
+ // At the end of the function, we return
+ // the norm of the residual:
+ return residual.l2_norm();
+ }
+
- // First thing to do after
- // refining the mesh, is to
- // setup the vectors,
- // matrices, etc., which is
- // done in the
- // <code>setup_system</code>
- // function.
- setup_system();
+ // @sect4{MinimalSurfaceProblem::determine_step_length}
- if (first_step)
- set_boundary_values ();
+ // As discussed in the introduction,
+ // Newton's method frequently does not
+ // converge if we always take full steps,
+ // i.e., compute $u^{n+1}=u^n+\delta
+ // u^n$. Rather, one needs a damping
+ // parameter (step length) $\alpha^n$ and
+ // set $u^{n+1}=u^n+\alpha^n\; delta
+ // u^n$. This function is the one called to
+ // compute $\alpha^n$.
+ //
+ // Here, we simply always return 0.1. This
+ // is of course a sub-optimal choice:
+ // ideally, what one wants is that the step
+ // size goes to one as we get closer to the
+ // solution, so that we get to enjoy the
+ // rapid quadratic convergence of Newton's
+ // method. We will discuss better
+ // strategies below in the results section.
+ template <int dim>
+ double MinimalSurfaceProblem<dim>::determine_step_length() const
+ {
+ return 0.1;
+ }
+
+
+
+ // @sect4{MinimalSurfaceProblem::run}
+
+ // In the run function, we build the first
+ // grid and then have the top-level logic
+ // for the Newton iteration. The function
+ // has two variables, one that indicates
+ // whether this is the first time we solve
+ // for a Newton update and one that
+ // indicates the refinement level of the
+ // mesh:
+ template <int dim>
+ void MinimalSurfaceProblem<dim>::run ()
+ {
+ unsigned int refinement = 0;
+ bool first_step = true;
+
+ // As described in the introduction, the
+ // domain is the unit disk around the
+ // origin, created in the same way as
+ // shown in step-6. The mesh is globally
+ // refined twice followed later on by
+ // several adaptive cycles:
+ GridGenerator::hyper_ball (triangulation);
+ static const HyperBallBoundary<dim> boundary;
+ triangulation.set_boundary (0, boundary);
+ triangulation.refine_global(2);
- // On every mesh there are
- // done five Newton steps, in
- // order to get a better
- // solution, before the mesh
- // gets too fine and the
- // computations take more
- // time.
- std::cout<<"initial residual:"<<compute_residual(0)<<std::endl;
+ // The Newton iteration starts
+ // next. During the first step we do not
+ // have information about the residual
+ // prior to this step and so we continue
+ // the Newton iteration until we have
+ // reached at least one iteration and
+ // until residual is less than $10^{-3}$.
+ //
+ // At the beginning of the loop, we do a
+ // bit of setup work. In the first go
+ // around, we compute the solution on the
+ // twice globally refined mesh after
+ // setting up the basic data
+ // structures. In all following mesh
+ // refinement loops, the mesh will be
+ // refined adaptively.
+ double previous_res = 0;
+ while (first_step || (previous_res>1e-3))
+ {
+ if (first_step == true)
+ {
+ std::cout << "******** Initial mesh "
+ << " ********"
+ << std::endl;
- for(unsigned int i=0; i<5;++i)
+ setup_system (true);
+ set_boundary_values ();
+ }
+ else
{
+ ++refinement;
+ std::cout << "******** Refined mesh " << refinement
+ << " ********"
+ << std::endl;
- // In every Newton step
- // the system matrix and
- // the right hand side
- // have to be computed.
+ refine_mesh();
+ }
+ // On every mesh we do exactly five
+ // Newton steps. We print the initial
+ // residual here and then start the
+ // iterations on this mesh.
+ //
+ // In every Newton step the system
+ // matrix and the right hand side
+ // have to be computed first, after
+ // which we store the norm of the
+ // right hand side as the residual to
+ // check against when deciding
+ // whether to stop the iterations. We
+ // then solve the linear system (the
+ // function also updates
+ // $u^{n+1}=u^n+\alpha^n\;\delta
+ // u^n$) and output the residual at
+ // the end of this Newton step:
+ std::cout << " Initial residual: "
+ << compute_residual(0)
+ << std::endl;
+
+ for (unsigned int inner_iteration=0; inner_iteration<5; ++inner_iteration)
+ {
assemble_system ();
previous_res = system_rhs.l2_norm();
solve ();
- first_step=false;
- std::cout<<"residual:"<<compute_residual(0)<<std::endl;
- }
- // The fifth solution, as
- // well as the Newton update,
- // on every mesh will be
- // written in a vtk-file, in
- // order to show the
- // convergence of the
- // solution.
-
- Assert (refinement < 100, ExcNotImplemented());
+ first_step = false;
+ std::cout << " Residual: "
+ << compute_residual(0)
+ << std::endl;
+ }
+ // Every fifth iteration, i.e., just
+ // before we refine the mesh again,
+ // we output the solution as well as
+ // the Newton update. This happens as
+ // in all programs before:
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
- data_out.add_data_vector (newton_update, "update");
data_out.add_data_vector (present_solution, "solution");
+ data_out.add_data_vector (newton_update, "update");
data_out.build_patches ();
- const std::string filename = "solution-" + Utilities::int_to_string (refinement, 2) + ".vtk";
+ const std::string filename = "solution-" +
+ Utilities::int_to_string (refinement, 2) +
+ ".vtk";
std::ofstream output (filename.c_str());
data_out.write_vtk (output);
// @sect4{The main function}
- // Finally the main function, this
+ // Finally the main function. This
// follows the scheme of all other
// main functions:
int main ()
projects / dealii-svn.git / commitdiff
