From 8d21cd48f30313ae0acc88e7e2f7dd8a934decb2 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 16 May 2012 12:47:00 +0000 Subject: [PATCH] Clean up a tiny bit. git-svn-id: https://svn.dealii.org/trunk@25511 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-15/step-15.cc | 934 ++++++++++++++-------------- 1 file changed, 472 insertions(+), 462 deletions(-) diff --git a/deal.II/examples/step-15/step-15.cc b/deal.II/examples/step-15/step-15.cc index c5bfa569b8..e0d0d63225 100644 --- a/deal.II/examples/step-15/step-15.cc +++ b/deal.II/examples/step-15/step-15.cc @@ -1,9 +1,9 @@ /* $Id$ */ -/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */ +/* Author: Sven Wetterauer, University of Heidelberg, 2012 */ /* $Id$ */ /* */ -/* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2010, 2011 by the deal.II authors */ +/* Copyright (C) 2012 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -12,31 +12,36 @@ // @sect3{Include files} - // The first few files have already - // been covered in previous examples - // and will thus not be further - // commented on. + // The first few files have already + // been covered in previous examples + // and will thus not be further + // commented on. #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 @@ -46,594 +51,599 @@ #include #include - // 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 to the new mesh. #include - // In this tutorial, we can't use the CG-method as a solver, as - // described above, but we use the minimal residual method, which - // is included with this file. - -#include + // In this tutorial, we can't use the CG-method as a solver, as + // described above, but we use the minimal residual method, which + // is included with this file. +#include +//#include - // As in previous programs: - -using namespace dealii; + // As in previous programs: +namespace Step15 +{ + using namespace dealii; - // @sect3{The Step15 class template} + // @sect3{The MinimalSurfaceProblem 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, 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. -template -class Step15 -{ - public: - Step15 (); - ~Step15 (); + template + class MinimalSurfaceProblem + { + public: + MinimalSurfaceProblem (); + ~MinimalSurfaceProblem (); - void run (); + void run (); - private: - void setup_system (); - void assemble_system (); - void solve (); - void refine_grid (); - void set_boundary_values (); - double compute_residual (const double alpha) const; - double determine_step_length() const; + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void set_boundary_values (); + double compute_residual (const double alpha) const; + double determine_step_length() const; - Triangulation triangulation; + Triangulation triangulation; - DoFHandler dof_handler; - FE_Q fe; + DoFHandler dof_handler; + FE_Q fe; - ConstraintMatrix hanging_node_constraints; + ConstraintMatrix hanging_node_constraints; - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; - Vector present_solution; - Vector newton_update; - Vector system_rhs; + Vector present_solution; + Vector newton_update; + Vector system_rhs; - unsigned int refinement; + 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; -}; + // 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} + // @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. + // 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. -template -class BoundaryValues : public Function -{ - public: - BoundaryValues () : Function() {} + template + class BoundaryValues : public Function + { + public: + BoundaryValues () : Function() {} - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; -template -double BoundaryValues::value (const Point &p, - const unsigned int /*component*/) const -{ - double return_value=sin(2*M_PI*(p[0]+p[1])); - return return_value; -} + template + double BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const + { + double return_value=sin(2*M_PI*(p[0]+p[1])); + return return_value; + } - // @sect3{The Step15 class implementation} + // @sect3{The MinimalSurfaceProblem class implementation} - // @sect4{Step15::Step15} + // @sect4{MinimalSurfaceProblem::MinimalSurfaceProblem} - // The constructor and destructor of the class are the same - // as in the first few tutorials. + // The constructor and destructor of the class are the same + // as in the first few tutorials. -template -Step15::Step15 () - : - dof_handler (triangulation), - fe (2) -{} + template + MinimalSurfaceProblem::MinimalSurfaceProblem () + : + dof_handler (triangulation), + fe (2) + {} - // -template -Step15::~Step15 () -{ - dof_handler.clear (); -} - - // @sect4{Step15::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. + // + template + MinimalSurfaceProblem::~MinimalSurfaceProblem () + { + dof_handler.clear (); + } -template -void Step15::setup_system () -{ + // @sect4{MinimalSurfaceProblem::setup_system} - // 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 refine_grid function, while the - // vector is transfered to the new mesh. + // 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. - if(first_step) + template + void MinimalSurfaceProblem::setup_system () { + + // 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 refine_grid function, while the + // vector is transfered to the new mesh. + + if (first_step) + { dof_handler.distribute_dofs (fe); present_solution.reinit (dof_handler.n_dofs()); for(unsigned int i=0; irefine_grid function - // after refining the mesh. + { + present_solution(i)=0; + } + // The constraint matrix, holding a list of the hanging nodes, + // will be setup in the refine_grid function + // after refining the mesh. hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, - hanging_node_constraints); + hanging_node_constraints); hanging_node_constraints.close (); + } + + + // The remaining parts of the function are the same as in step 6. + + newton_update.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, c_sparsity); + + hanging_node_constraints.condense (c_sparsity); + + sparsity_pattern.copy_from(c_sparsity); + system_matrix.reinit (sparsity_pattern); } + // @sect4{MinimalSurfaceProblem::assemble_system} - // The remaining parts of the function are the same as in step 6. + // 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. - newton_update.reinit (dof_handler.n_dofs()); - system_rhs.reinit (dof_handler.n_dofs()); + template + void MinimalSurfaceProblem::assemble_system () + { + const QGauss quadrature_formula(3); - CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, c_sparsity); + system_matrix = 0; + system_rhs = 0; - hanging_node_constraints.condense (c_sparsity); + FEValues fe_values (fe, quadrature_formula, + update_gradients | + update_quadrature_points | update_JxW_values); - sparsity_pattern.copy_from(c_sparsity); - system_matrix.reinit (sparsity_pattern); -} + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); - // @sect4{Step15::assemble_system} + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); - // 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. + std::vector local_dof_indices (dofs_per_cell); -template -void Step15::assemble_system () -{ - const QGauss quadrature_formula(3); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; - system_matrix = 0; - system_rhs = 0; + fe_values.reinit (cell); - FEValues 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(); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); + // 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 local_dof_indices (dofs_per_cell); + std::vector > gradients(n_q_points); + fe_values.get_function_gradients(present_solution, gradients); - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - cell_matrix = 0; - 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 > 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 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 - * (fe_values.shape_grad(j, q_point) - * gradients[q_point]) - * 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->get_dof_indices (local_dof_indices); - for (unsigned int i=0; i boundary_values; + + cell->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction(), - boundary_values); - - MatrixTools::apply_boundary_values (boundary_values, - system_matrix, - newton_update, - system_rhs); -} + 0, + ZeroFunction(), + boundary_values); + + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + newton_update, + system_rhs); + } -template -double Step15::compute_residual (const double alpha) const -{ - const QGauss quadrature_formula(3); + template + double MinimalSurfaceProblem::compute_residual (const double alpha) const + { + const QGauss quadrature_formula(3); - Vector residual (dof_handler.n_dofs()); + Vector residual (dof_handler.n_dofs()); - Vector linearization_point (dof_handler.n_dofs()); - linearization_point = present_solution; - linearization_point.add (alpha, newton_update); + Vector linearization_point (dof_handler.n_dofs()); + linearization_point = present_solution; + linearization_point.add (alpha, newton_update); - FEValues fe_values (fe, quadrature_formula, - update_gradients | - update_quadrature_points | update_JxW_values); + FEValues 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(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); - Vector cell_rhs (dofs_per_cell); + Vector cell_rhs (dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - cell_rhs = 0; + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_rhs = 0; - fe_values.reinit (cell); + fe_values.reinit (cell); - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { + 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: + // 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 > gradients(n_q_points); - fe_values.get_function_gradients(linearization_point, gradients); + std::vector > 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. + // 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]); + 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. + // 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)); - } - } + 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 boundary_values; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction(), - boundary_values); - for (std::map::const_iterator p = boundary_values.begin(); - p != boundary_values.end(); ++p) - residual(p->first) = 0; - - return residual.l2_norm(); -} + cell->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + boundary_values); + for (std::map::const_iterator p = boundary_values.begin(); + p != boundary_values.end(); ++p) + residual(p->first) = 0; + + return residual.l2_norm(); + } - // 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. + // @sect4{MinimalSurfaceProblem::solve} -template -void Step15::solve () -{ - SolverControl solver_control (1000, system_rhs.l2_norm()*1e-6); - SolverMinRes<> solver (solver_control); + // 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. + + template + void MinimalSurfaceProblem::solve () + { + SolverControl solver_control (system_rhs.size(), + system_rhs.l2_norm()*1e-6); + SolverCG<> solver (solver_control); - PreconditionSSOR<> preconditioner; - preconditioner.initialize(system_matrix, 1.2); + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); - solver.solve (system_matrix, newton_update, system_rhs, - preconditioner); + solver.solve (system_matrix, newton_update, system_rhs, + preconditioner); - hanging_node_constraints.distribute (newton_update); + 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); -} + // 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 -double Step15::determine_step_length() const -{ - return 0.1; -} - // @sect4{Step15::refine_grid} + template + double MinimalSurfaceProblem::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. + // 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. -template -void Step15::refine_grid () -{ - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - KellyErrorEstimator::estimate (dof_handler, - QGauss(3), - typename FunctionMap::type(), - present_solution, - estimated_error_per_cell); - - GridRefinement::refine_and_coarsen_fixed_number (triangulation, - 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 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 enforce rules like the 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. - - 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: - - SolutionTransfer 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: - - Vector tmp(dof_handler.n_dofs()); - 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: - hanging_node_constraints.clear(); - - // After doing so, the hanging nodes of the new mesh can be - // enlisted in the matrix, like before. Calling the - // setup_system function in the run - // 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); - hanging_node_constraints.close(); - hanging_node_constraints.distribute(present_solution); -} + template + void MinimalSurfaceProblem::refine_grid () + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::estimate (dof_handler, + QGauss(3), + typename FunctionMap::type(), + present_solution, + estimated_error_per_cell); + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + 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 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 enforce rules like the 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. + + 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: + + SolutionTransfer 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: + + Vector tmp(dof_handler.n_dofs()); + 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: + hanging_node_constraints.clear(); + + // After doing so, the hanging nodes of the new mesh can be + // enlisted in the matrix, like before. Calling the + // setup_system function in the run + // 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); + hanging_node_constraints.close(); + hanging_node_constraints.distribute(present_solution); + } -template -void Step15::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: + template + void MinimalSurfaceProblem::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 boundary_values2; VectorTools::interpolate_boundary_values(dof_handler, 0, - BoundaryValues(), boundary_values2); + BoundaryValues(), boundary_values2); for (std::map::const_iterator p = - boundary_values2.begin(); p != boundary_values2.end(); ++p) + boundary_values2.begin(); p != boundary_values2.end(); ++p) present_solution(p->first) = p->second; -} - // @sect4{Step15::run} + } + // @sect4{MinimalSurfaceProblem::run} - // In the run function, the first grid is build. Also in this - // function, the Newton iteration is implemented. + // In the run function, the first grid is build. Also in this + // function, the Newton iteration is implemented. -template -void Step15::run () -{ + template + void MinimalSurfaceProblem::run () + { - // The integer refinement counts the mesh refinements. Obviously - // starting the program, it should be zero. - refinement=0; - first_step=true; + // 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. + // 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. - GridGenerator::hyper_ball (triangulation); - static const HyperBallBoundary boundary; - triangulation.set_boundary (0, boundary); - triangulation.refine_global(2); + GridGenerator::hyper_ball (triangulation); + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + triangulation.refine_global(2); - // 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}$. + // 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}$. - double previous_res = 0; - while(first_step || (previous_res>1e-3)) - { + double previous_res = 0; + while(first_step || (previous_res>1e-3)) + { - // 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) - { - refine_grid(); + // 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) + { + refine_grid(); - std::cout<<"********mesh-refinement:"<setup_system - // function. + // First thing to do after refining the mesh, is to setup the vectors, + // matrices, etc., which is done in the setup_system + // function. - setup_system(); + setup_system(); - if (first_step) - set_boundary_values (); + if (first_step) + set_boundary_values (); - // 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:"< data_out; + DataOut 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.build_patches (); - const std::string filename = "solution-" + Utilities::int_to_string (refinement, 2) + ".vtk"; - std::ofstream output (filename.c_str()); - data_out.write_vtk (output); + 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.build_patches (); + 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} + // @ sect4{The main function} - // Finally the main function, this follows the scheme of all other main - // functions: + // Finally the main function, this follows the scheme of all other main + // functions: int main () { - try { + using namespace dealii; + using namespace Step15; + deallog.depth_console (0); - Step15<2> laplace_problem_2d; + MinimalSurfaceProblem<2> laplace_problem_2d; laplace_problem_2d.run (); } catch (std::exception &exc) -- 2.39.5