From 54f039681b6c8d92c289569cd33255bc1f5bbd6d Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 9 Feb 2012 13:32:20 +0000 Subject: [PATCH] Minor reshuffling. git-svn-id: https://svn.dealii.org/trunk@25021 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-41/step-41.cc | 1020 ++++++++++++++------------- 1 file changed, 525 insertions(+), 495 deletions(-) diff --git a/deal.II/examples/step-41/step-41.cc b/deal.II/examples/step-41/step-41.cc index 131b883ae4..911404cf08 100644 --- a/deal.II/examples/step-41/step-41.cc +++ b/deal.II/examples/step-41/step-41.cc @@ -1,9 +1,8 @@ -/* $Id: step-4.cc 24093 2011-08-16 13:58:12Z bangerth $ */ -/* Author: Wolfgang Bangerth, University of Heidelberg, 1999 */ +/* Author: Joerg Frohne Wolfgang Bangerth, Texas A&M University, 2011, 2012 */ /* $Id: step-4.cc 24093 2011-08-16 13:58:12Z bangerth $ */ /* */ -/* Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors */ +/* Copyright (C) 2011, 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,544 +11,545 @@ // @sect3{Include files} -#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 - -using namespace dealii; - - // @sect3{The Step41 class template} - - // This class supply all function and variables - // to an obstacle problem. The projection_active_set - // function and the ConstaintMatrix are important - // for the handling of the active set as we see - // later. - -template -class Step41 -{ - public: - Step41 (); - void run (); - - private: - void make_grid (); - void setup_system(); - void assemble_system (); - void assemble_mass_matrix (); - void projection_active_set (); - void solve (); - void output_results (const std::string& title) const; - - Triangulation triangulation; - FE_Q fe; - DoFHandler dof_handler; - unsigned int n_refinements; - - ConstraintMatrix constraints; - - SparsityPattern sparsity_pattern; - TrilinosWrappers::SparseMatrix system_matrix; - TrilinosWrappers::SparseMatrix system_matrix_complete; - TrilinosWrappers::SparseMatrix mass_matrix; - - TrilinosWrappers::Vector solution; - TrilinosWrappers::Vector tmp_solution; - TrilinosWrappers::Vector system_rhs; - TrilinosWrappers::Vector system_rhs_complete; - TrilinosWrappers::Vector resid_vector; - TrilinosWrappers::Vector active_set; - TrilinosWrappers::Vector diag_mass_matrix_vector; -}; - - - // @sect3{Right hand side and boundary values} - -template -class RightHandSide : public Function -{ - public: - RightHandSide () : Function() {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - -template -class BoundaryValues : public Function -{ - public: - BoundaryValues () : Function() {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - -template -class Obstacle : public Function -{ - public: - Obstacle () : Function() {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; +#include +#include +#include +#include +#include +#include +#include +#include +#include - // For this example, we choose as right hand - // side function a constant force density - // like the gravitation attraction. -template -double RightHandSide::value (const Point &p, - const unsigned int /*component*/) const -{ - double return_value = -10; +#include +#include - return return_value; -} +#include +#include +#include - // As boundary values, we choose the zero. -template -double BoundaryValues::value (const Point &p, - const unsigned int /*component*/) const +namespace Step41 { - double return_value = 0; + using namespace dealii; + + // @sect3{The Step41 class template} + + // This class supply all function and variables + // to an obstacle problem. The projection_active_set + // function and the ConstaintMatrix are important + // for the handling of the active set as we see + // later. + + template + class ObstacleProblem + { + public: + ObstacleProblem (); + void run (); + + private: + void make_grid (); + void setup_system(); + void assemble_system (); + void assemble_mass_matrix (); + void projection_active_set (); + void solve (); + void output_results (const std::string& title) const; + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + unsigned int n_refinements; + + ConstraintMatrix constraints; + + SparsityPattern sparsity_pattern; + TrilinosWrappers::SparseMatrix system_matrix; + TrilinosWrappers::SparseMatrix system_matrix_complete; + TrilinosWrappers::SparseMatrix mass_matrix; + + TrilinosWrappers::Vector solution; + TrilinosWrappers::Vector tmp_solution; + TrilinosWrappers::Vector system_rhs; + TrilinosWrappers::Vector system_rhs_complete; + TrilinosWrappers::Vector resid_vector; + TrilinosWrappers::Vector active_set; + TrilinosWrappers::Vector diag_mass_matrix_vector; + }; + + + // @sect3{Right hand side and boundary values} + + template + class RightHandSide : public Function + { + public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; + + template + class BoundaryValues : public Function + { + public: + BoundaryValues () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; + + template + class Obstacle : public Function + { + public: + Obstacle () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; + + + + // For this example, we choose as right hand + // side function a constant force density + // like the gravitation attraction. + template + double RightHandSide::value (const Point &p, + const unsigned int component) const + { + Assert (component == 0, ExcNotImplemented()); + + return -10; + } + + + // As boundary values, we choose the zero. + template + double BoundaryValues::value (const Point &p, + const unsigned int component) const + { + Assert (component == 0, ExcNotImplemented()); + + return 0; + } + + + // The obstacle function describes a cascaded + // barrier. So if the gravitation attraction + // pulls the membrane down it blows over the + // steps. + template + double Obstacle::value (const Point &p, + const unsigned int component) const + { + Assert (component == 0, ExcNotImplemented()); + + if (p (0) < -0.5) + return -0.2; + else if (p (0) >= -0.5 && p (0) < 0.0) + return -0.4; + else if (p (0) >= 0.0 && p (0) < 0.5) + return -0.6; + else + return -0.8; + } - return return_value; -} - // The obstacle function describes a cascaded - // barrier. So if the gravitation attraction - // pulls the membrane down it blows over the - // steps. -template -double Obstacle::value (const Point &p, - const unsigned int /*component*/) const -{ - double return_value = 0; - - if (p (0) < -0.5) - return_value = -0.2; - else if (p (0) >= -0.5 && p (0) < 0.0) - return_value = -0.4; - else if (p (0) >= 0.0 && p (0) < 0.5) - return_value = -0.6; - else - return_value = -0.8; - - return return_value; -} + // @sect3{Implementation of the ObstacleProblem class} + + + // @sect4{ObstacleProblem::ObstacleProblem} + + template + ObstacleProblem::ObstacleProblem () + : + fe (1), + dof_handler (triangulation) + {} + // @sect4{ObstacleProblem::make_grid} - // @sect3{Implementation of the Step41 class} - + // We solve our obstacle problem on the square + // $[-1,1]\times [-1,1]$ in 2D. + template + void ObstacleProblem::make_grid () + { + GridGenerator::hyper_cube (triangulation, -1, 1); + n_refinements = 8; + triangulation.refine_global (n_refinements); + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl + << " Total number of cells: " + << triangulation.n_cells() + << std::endl; + } - // @sect4{Step41::Step41} + // @sect4{ObstacleProblem::setup_system} -template -Step41::Step41 () - : - fe (1), - dof_handler (triangulation) -{} + template + void ObstacleProblem::setup_system () + { + dof_handler.distribute_dofs (fe); + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; - // @sect4{Step41::make_grid} + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, c_sparsity, constraints, false); + sparsity_pattern.copy_from(c_sparsity); - // We solve our obstacle problem on the square - // $[-1,1]\times [-1,1]$ in 2D. -template -void Step41::make_grid () -{ - GridGenerator::hyper_cube (triangulation, -1, 1); - n_refinements = 6; - triangulation.refine_global (n_refinements); - - std::cout << " Number of active cells: " - << triangulation.n_active_cells() - << std::endl - << " Total number of cells: " - << triangulation.n_cells() - << std::endl; -} + system_matrix.reinit (sparsity_pattern); + system_matrix_complete.reinit (sparsity_pattern); + mass_matrix.reinit (sparsity_pattern); - // @sect4{Step41::setup_system} + solution.reinit (dof_handler.n_dofs()); + tmp_solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + system_rhs_complete.reinit (dof_handler.n_dofs()); + resid_vector.reinit (dof_handler.n_dofs()); + active_set.reinit (dof_handler.n_dofs()); + diag_mass_matrix_vector.reinit (dof_handler.n_dofs()); + } -template -void Step41::setup_system () -{ - dof_handler.distribute_dofs (fe); - - std::cout << " Number of degrees of freedom: " - << dof_handler.n_dofs() - << std::endl; - - CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, c_sparsity, constraints, false); - sparsity_pattern.copy_from(c_sparsity); - - system_matrix.reinit (sparsity_pattern); - system_matrix_complete.reinit (sparsity_pattern); - mass_matrix.reinit (sparsity_pattern); - - solution.reinit (dof_handler.n_dofs()); - tmp_solution.reinit (dof_handler.n_dofs()); - system_rhs.reinit (dof_handler.n_dofs()); - system_rhs_complete.reinit (dof_handler.n_dofs()); - resid_vector.reinit (dof_handler.n_dofs()); - active_set.reinit (dof_handler.n_dofs()); - diag_mass_matrix_vector.reinit (dof_handler.n_dofs()); -} + // @sect4{ObstacleProblem::assemble_system} - // @sect4{Step41::assemble_system} + // At once with assembling the system matrix and + // right-hand-side we apply the constraints + // to our system. The constraint consists not + // only of the zero Dirichlet boundary values, + // in addition they contain the obstacle values. + // The projection_active_set function are used + // to fill the ConstraintMatrix. + template + void ObstacleProblem::assemble_system () + { + QGauss quadrature_formula(2); - // At once with assembling the system matrix and - // right-hand-side we apply the constraints - // to our system. The constraint consists not - // only of the zero Dirichlet boundary values, - // in addition they contain the obstacle values. - // The projection_active_set function are used - // to fill the ConstraintMatrix. -template -void Step41::assemble_system () -{ - QGauss quadrature_formula(2); + const RightHandSide right_hand_side; - const RightHandSide right_hand_side; + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); - FEValues fe_values (fe, quadrature_formula, - update_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 cell_matrix (dofs_per_cell, dofs_per_cell); + TrilinosWrappers::Vector cell_rhs (dofs_per_cell); - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - TrilinosWrappers::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(); - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - cell_matrix = 0; - cell_rhs = 0; - - for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - - // This function apply the constraints - // to the system matrix and system rhs. - // The true parameter is set to make sure - // that the system rhs contains correct - // values in the rows with inhomogeneity - // constraints. - constraints.distribute_local_to_global (cell_matrix, cell_rhs, - local_dof_indices, - system_matrix, system_rhs, true); - } -} + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_matrix = 0; + cell_rhs = 0; + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + + // This function apply the constraints + // to the system matrix and system rhs. + // The true parameter is set to make sure + // that the system rhs contains correct + // values in the rows with inhomogeneity + // constraints. + constraints.distribute_local_to_global (cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, system_rhs, true); + } + } -template -void Step41::assemble_mass_matrix () -{ - QTrapez quadrature_formula; + template + void ObstacleProblem::assemble_mass_matrix () + { + QTrapez quadrature_formula; - FEValues fe_values (fe, quadrature_formula, - update_values | update_quadrature_points | update_JxW_values); + FEValues fe_values (fe, quadrature_formula, + update_values | 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 cell_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix cell_matrix (dofs_per_cell, 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) - { - fe_values.reinit (cell); - cell_matrix = 0; - - for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - - // This function apply the constraints - // to the system matrix and system rhs. - // The true parameter is set to make sure - // that the system rhs contains correct - // values in the rows with inhomogeneity - // constraints. - constraints.distribute_local_to_global (cell_matrix, - local_dof_indices, - mass_matrix); - } -} + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); - // @sect4{Step41::projection_active_set} - - // Updating of the active set which means to - // set a inhomogeneity constraint in the - // ConstraintMatrix. At the same time we set - // the solution to the correct value - the obstacle value. - // To control the active set we use the vector - // active_set which contains a zero in a component - // that is not in the active set and elsewise a - // one. With the output file you can visualize it. -template -void Step41::projection_active_set () -{ - const Obstacle obstacle; - std::vector vertex_touched (triangulation.n_vertices(), - false); - unsigned int counter_contact_constraints = 0; - - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - constraints.clear(); - - // to find and supply the constraints for the - // obstacle condition - active_set = 0.0; - const double c = 100.0; - for (; cell!=endc; ++cell) - for (unsigned int v=0; v::vertices_per_cell; ++v) + for (; cell!=endc; ++cell) { - unsigned int index_x = cell->vertex_dof_index (v,0); + fe_values.reinit (cell); + cell_matrix = 0; - // the local row where - Point point (cell->vertex (v)[0], cell->vertex (v)[1]); - double obstacle_value = obstacle.value (point); - double solution_index_x = solution (index_x); - - // To decide which dof belongs to the - // active-set. For that we scale the - // residual-vector with the cell-size and - // the diag-entry of the mass-matrix. + for (unsigned int q_point=0; q_point 0) + cell->get_dof_indices (local_dof_indices); + + // This function apply the constraints + // to the system matrix and system rhs. + // The true parameter is set to make sure + // that the system rhs contains correct + // values in the rows with inhomogeneity + // constraints. + constraints.distribute_local_to_global (cell_matrix, + local_dof_indices, + mass_matrix); + } + } + + // @sect4{ObstacleProblem::projection_active_set} + + // Updating of the active set which means to + // set a inhomogeneity constraint in the + // ConstraintMatrix. At the same time we set + // the solution to the correct value - the obstacle value. + // To control the active set we use the vector + // active_set which contains a zero in a component + // that is not in the active set and elsewise a + // one. With the output file you can visualize it. + template + void ObstacleProblem::projection_active_set () + { + const Obstacle obstacle; + std::vector vertex_touched (triangulation.n_vertices(), + false); + unsigned int counter_contact_constraints = 0; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + constraints.clear(); + + // to find and supply the constraints for the + // obstacle condition + active_set = 0.0; + const double c = 100.0; + for (; cell!=endc; ++cell) + for (unsigned int v=0; v::vertices_per_cell; ++v) { - constraints.add_line (index_x); - constraints.set_inhomogeneity (index_x, obstacle_value); - solution (index_x) = obstacle_value; - active_set (index_x) = 1.0; - - if (vertex_touched[cell->vertex_index(v)] == false) - { - vertex_touched[cell->vertex_index(v)] = true; - counter_contact_constraints += 1; - } + unsigned int index_x = cell->vertex_dof_index (v,0); + + // the local row where + Point point (cell->vertex (v)[0], cell->vertex (v)[1]); + double obstacle_value = obstacle.value (point); + double solution_index_x = solution (index_x); + + // To decide which dof belongs to the + // active-set. For that we scale the + // residual-vector with the cell-size and + // the diag-entry of the mass-matrix. + + // TODO: I have to check the condition + if (resid_vector (index_x) + + diag_mass_matrix_vector (index_x)*c*(obstacle_value - solution_index_x) > 0) + { + constraints.add_line (index_x); + constraints.set_inhomogeneity (index_x, obstacle_value); + solution (index_x) = obstacle_value; + active_set (index_x) = 1.0; + + if (vertex_touched[cell->vertex_index(v)] == false) + { + vertex_touched[cell->vertex_index(v)] = true; + counter_contact_constraints += 1; + } + } } + std::cout<< "Number of Contact-Constaints: " << counter_contact_constraints <(), + constraints); + constraints.close (); + } + + // @sect4{ObstacleProblem::solve} + + template + void ObstacleProblem::solve () + { + ReductionControl reduction_control (100, 1e-12, 1e-3); + SolverCG solver (reduction_control); + TrilinosWrappers::PreconditionAMG precondition; + precondition.initialize (system_matrix); + + solver.solve (system_matrix, solution, system_rhs, precondition); + + std::cout << "Initial error: " << reduction_control.initial_value() < + void ObstacleProblem::output_results (const std::string& title) const + { + DataOut data_out; + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (tmp_solution, "Displacement"); + data_out.add_data_vector (resid_vector, "Residual"); + data_out.add_data_vector (active_set, "ActiveSet"); + + data_out.build_patches (); + + std::ofstream output_vtk ((title + ".vtk").c_str ()); + data_out.write_vtk (output_vtk); + } + + + + // @sect4{ObstacleProblem::run} + + // This is the function which has the + // top-level control over everything. + // Here the active set method is implemented. + + // TODO: I have to compare it with the algorithm + // in the Wohlmuth-paper + template + void ObstacleProblem::run () + { + std::cout << "Solving problem in " << dim << " space dimensions." << std::endl; + + make_grid(); + setup_system (); + + constraints.clear (); + VectorTools::interpolate_boundary_values (dof_handler, + 0, + BoundaryValues(), + constraints); + constraints.close (); + ConstraintMatrix constraints_complete (constraints); + assemble_system (); + solve (); + + // to save the system_matrix and the + // rhs to compute the residual in every + // step of the active-set-iteration + system_matrix_complete.copy_from (system_matrix); + system_rhs_complete = system_rhs; + + // to compute the factor which is used + // to scale the residual. You can consider + // this diagonal matrix as the discretization + // of a lagrange multiplier for the + // contact force + assemble_mass_matrix (); + for (unsigned int j=0; j(), - constraints); - constraints.close (); -} - - // @sect4{Step41::solve} - -template -void Step41::solve () -{ - ReductionControl reduction_control (100, 1e-12, 1e-3); - SolverCG solver (reduction_control); - TrilinosWrappers::PreconditionAMG precondition; - precondition.initialize (system_matrix); - - solver.solve (system_matrix, solution, system_rhs, precondition); - - std::cout << "Initial error: " << reduction_control.initial_value() < -void Step41::output_results (const std::string& title) const -{ - DataOut data_out; - - data_out.attach_dof_handler (dof_handler); - // data_out.add_data_vector (tmp_solution, "Displacement"); - // data_out.add_data_vector (resid_vector, "Residual"); - data_out.add_data_vector (active_set, "ActiveSet"); - - data_out.build_patches (); - - std::ofstream output_vtk ((title + ".vtk").c_str ()); - data_out.write_gnuplot (output_vtk); -} - - - - // @sect4{Step41::run} - - // This is the function which has the - // top-level control over everything. - // Here the active set method is implemented. - - // TODO: I have to compare it with the algorithm - // in the Wohlmuth-paper -template -void Step41::run () -{ - std::cout << "Solving problem in " << dim << " space dimensions." << std::endl; - - make_grid(); - setup_system (); - - constraints.clear (); - VectorTools::interpolate_boundary_values (dof_handler, - 0, - BoundaryValues(), - constraints); - constraints.close (); - ConstraintMatrix constraints_complete (constraints); - assemble_system (); - solve (); - - // to save the system_matrix and the - // rhs to compute the residual in every - // step of the active-set-iteration - system_matrix_complete.copy_from (system_matrix); - system_rhs_complete = system_rhs; - - // to compute the factor which is used - // to scale the residual. You can consider - // this diagonal matrix as the discretization - // of a lagrange multiplier for the - // contact force - assemble_mass_matrix (); - for (unsigned int j=0; j laplace_problem_2d; - laplace_problem_2d.run (); - return 0; } -- 2.39.5