From: sfloresm Date: Mon, 9 Nov 2015 13:53:16 +0000 (-0300) Subject: Solve p-laplacian elliptic problems for large p X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a4a37bda88d77898c852c172d4f9f7a68a8da2ac;p=code-gallery.git Solve p-laplacian elliptic problems for large p --- diff --git a/Elastoplastic Torsion/CMakeLists.txt b/Elastoplastic Torsion/CMakeLists.txt deleted file mode 100644 index 82fb113..0000000 --- a/Elastoplastic Torsion/CMakeLists.txt +++ /dev/null @@ -1,40 +0,0 @@ -## -# CMake script for the Elastoplastic Torsion problem -# at the deal.ii code-gallery -## - -# Set the name of the project and target: -SET(TARGET "ElastoplasticTorsion") - -# Declare all source files the target consists of. Here, this is only -# the one step-X.cc file, but as you expand your project you may wish -# to add other source files as well. If your project becomes much larger, -# you may want to either replace the following statement by something like -# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") -# FILE(GLOB_RECURSE TARGET_INC "include/*.h") -# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) -# or switch altogether to the large project CMakeLists.txt file discussed -# in the "CMake in user projects" page accessible from the "User info" -# page of the documentation. -SET(TARGET_SRC - ${TARGET}.cc - ) - -# Usually, you will not need to modify anything beyond this point... - -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) - -FIND_PACKAGE(deal.II 8.0 QUIET - HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} - ) -IF(NOT ${deal.II_FOUND}) - MESSAGE(FATAL_ERROR "\n" - "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" - "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" - "or set an environment variable \"DEAL_II_DIR\" that contains this path." - ) -ENDIF() - -DEAL_II_INITIALIZE_CACHED_VARIABLES() -PROJECT(${TARGET}) -DEAL_II_INVOKE_AUTOPILOT() diff --git a/Elastoplastic Torsion/ElastoplasticTorsion.cc b/Elastoplastic Torsion/ElastoplasticTorsion.cc deleted file mode 100644 index 1cec33a..0000000 --- a/Elastoplastic Torsion/ElastoplasticTorsion.cc +++ /dev/null @@ -1,1343 +0,0 @@ -/* --------------------------------------------------------------------- - * - * Copyright (C) 2010 - 2015 by the deal.II authors - * and Salvador Flores. - * - * - * - * This is free software; you can use it, redistribute - * it, and/or modify it under the terms of the GNU Lesser General - * Public License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * The full text of the license can be found in the file LICENSE at - * the top level of the deal.II distribution. - * - * --------------------------------------------------------------------- - * - * Author: Salvador Flores, - * Center for Mathematical Modelling, - * Universidad de Chile, 2015. - */ - - - /* - This piece of software solves the elliptic p-laplacian - boundary-value problems: - - Min {∫ 1/2 W(|Du|²)+ 1/p |Du|^p -fu : u=g on ∂S } (1) - u - - for large values of p, which approximates (see Alvarez & Flores 2015) - - Min {∫ 1/2 W(|Du|²) -fu : |Du|<1 a.s. on S, u=g on ∂S } - u - - By default W(t)=t and S=unit disk. - - Large portions of this code are borrowed from the deal.ii tutorials - - step-15, step-29. - - For further details see the technical report available at - the documentation andmof at http://www.dim.uchile.cl/~sflores. - - */ - -// 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 - -#include -#include -#include -#include -#include - -#include -#include -#include - -#include - -// Open a namespace for this program and import everything from the -// dealii namespace into it. -namespace nsp -{ - using namespace dealii; - -// ********************************************************// - class ParameterReader : public Subscriptor -{ -public: -ParameterReader(ParameterHandler &); -void read_parameters(const std::string); -private: -void declare_parameters(); -ParameterHandler &prm; -}; -// Constructor -ParameterReader::ParameterReader(ParameterHandler ¶mhandler): -prm(paramhandler) -{} - - void ParameterReader::declare_parameters() -{ - -prm.enter_subsection ("Global Parameters"); -{ - prm.declare_entry("p", "100",Patterns::Double(2.1), - "Penalization parameter"); - prm.declare_entry("known_solution", "true",Patterns::Bool(), - "Whether the exact solution is known"); -} -prm.leave_subsection (); - -prm.enter_subsection ("Mesh & Refinement Parameters"); -{ - prm.declare_entry("Code for the domain", "0",Patterns::Integer(0,2), - "Number identifying the domain in which we solve the problem"); - prm.declare_entry("No of initial refinements", "4",Patterns::Integer(0), - "Number of global mesh refinement steps applied to initial coarse grid"); - prm.declare_entry("No of adaptive refinements", "8",Patterns::Integer(0), - "Number of global adaptive mesh refinements"); - prm.declare_entry("top_fraction_of_cells", "0.25",Patterns::Double(0), - "refinement threshold"); - prm.declare_entry("bottom_fraction_of_cells", "0.05",Patterns::Double(0), - "coarsening threshold"); -} -prm.leave_subsection (); - - -prm.enter_subsection ("Algorithm Parameters"); -{ - prm.declare_entry("Descent_direction", "0",Patterns::Integer(0,1), - "0: Preconditioned descent, 1: Newton Method"); - prm.declare_entry("init_p", "10",Patterns::Double(2), - "Initial p"); - prm.declare_entry("delta_p", "50",Patterns::Double(0), - "increase of p"); - prm.declare_entry("Max_CG_it", "1500",Patterns::Integer(1), - "Maximum Number of CG iterations"); - prm.declare_entry("CG_tol", "1e-10",Patterns::Double(0), - "Tolerance for CG iterations"); - prm.declare_entry("max_LS_it", "45",Patterns::Integer(1), - "Maximum Number of LS iterations"); - prm.declare_entry("line_search_tolerence", "1e-6",Patterns::Double(0), - "line search tolerance constant (c1 in Nocedal-Wright)"); - prm.declare_entry("init_step_length", "1e-2",Patterns::Double(0), - "initial step length in line-search"); - prm.declare_entry("Max_inner", "800",Patterns::Integer(1), - "Maximum Number of inner iterations"); - prm.declare_entry("eps", "1.0e-8",Patterns::Double(0), - "Threshold on norm of the derivative to declare optimality achieved"); - prm.declare_entry("hi_eps", "1.0e-9",Patterns::Double(0), - "Threshold on norm of the derivative to declare optimality achieved in highly refined mesh"); - prm.declare_entry("hi_th", "8",Patterns::Integer(0), - "Number of adaptive refinement before change convergence threshold"); -} -prm.leave_subsection (); - -} -void ParameterReader::read_parameters (const std::string parameter_file) -{ -declare_parameters(); -prm.read_input (parameter_file); -} - -// ******************************************************************************************// -// The solution of the elastoplastic torsion problem on the unit disk with rhs=4. - -template -class Solution : public Function -{ -public: -Solution () : Function() {} -virtual double value (const Point &pto, const unsigned int component = 0) const; -virtual Tensor<1,dim> gradient (const Point &pto, const unsigned int component = 0) const; -}; - -template -double Solution::value (const Point &pto,const unsigned int) const -{ -double r=sqrt(pto.square()); - if (r<0.5) - return -1.0*std::pow(r,2.0)+0.75; - else - return 1.0-r; -} - - - -template -Tensor<1,dim> Solution::gradient (const Point &pto,const unsigned int) const -{ -double r=sqrt(pto.square()); - if (r<0.5) - return -2.0*pto; - else - return -1.0*pto/r; -} - - - - -// ****************************************************************************************** // -/* Compute the Lagrange multiplier (as a derived quantity) */ - - -template -class ComputeMultiplier : public DataPostprocessor -{ - private: - double p; - public: - ComputeMultiplier (double pe); - - virtual - void compute_derived_quantities_scalar ( - const std::vector< double > & , - const std::vector< Tensor< 1, dim > > & , - const std::vector< Tensor< 2, dim > > & , - const std::vector< Point< dim > > & , - const std::vector< Point< dim > > & , - std::vector< Vector< double > > & - ) const; - - virtual std::vector get_names () const; - - virtual - std::vector - get_data_component_interpretation () const; - virtual UpdateFlags get_needed_update_flags () const; -}; - - - template - ComputeMultiplier::ComputeMultiplier (double pe): p(pe) - {} - - -template -void ComputeMultiplier::compute_derived_quantities_scalar( - const std::vector< double > & /*uh*/, - const std::vector< Tensor< 1, dim > > & duh, - const std::vector< Tensor< 2, dim > > & /*dduh*/, - const std::vector< Point< dim > > & /* normals*/, - const std::vector< Point< dim > > & /*evaluation_points*/, - std::vector< Vector< double > > & computed_quantities ) const -{ - const unsigned int n_quadrature_points = duh.size(); - - for (unsigned int q=0; q -std::vector -ComputeMultiplier::get_names() const -{ - std::vector solution_names; -solution_names.push_back ("Gradient norm"); -solution_names.push_back ("Lagrange multiplier"); - return solution_names; -} - - -template -UpdateFlags -ComputeMultiplier::get_needed_update_flags () const -{ - return update_gradients; -} - - - - template - std::vector -ComputeMultiplier:: get_data_component_interpretation () const - { - std::vector - interpretation; - // norm of the gradient - interpretation.push_back (DataComponentInterpretation::component_is_scalar); - // Lagrange multiplier - interpretation.push_back (DataComponentInterpretation::component_is_scalar); - return interpretation; -} - - - - - -// *************************************************************************************** // - template - class ElastoplasticTorsion - { - public: - ElastoplasticTorsion (ParameterHandler &); - ~ElastoplasticTorsion (); - void run (); - - private: - void setup_system (const bool initial_step); - void assemble_system (); - bool solve (const int inner_it); - void init_mesh (); - void refine_mesh (); - void set_boundary_values (); - double phi (const double alpha) const; - bool checkWolfe(double & alpha, double & phi_alpha) const; - bool determine_step_length (const int inner_it); - void print_it_message (const int counter, bool ks); - void output_results (unsigned int refinement) const; - void format_convergence_tables(); - void process_solution (const unsigned int cycle); - void process_multiplier (const unsigned int cycle,const int iter,double time); - double dual_error () const; - double dual_infty_error () const; - double W (double Du2) const; - double Wp (double Du2) const; - double G (double Du2) const; - - - - ParameterHandler &prm; - Triangulation triangulation; - DoFHandler dof_handler; - ConstraintMatrix hanging_node_constraints; - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - ConvergenceTable convergence_table; - ConvergenceTable dual_convergence_table; - Vector present_solution; - Vector newton_update; - Vector system_rhs; - Vector grad_norm; - Vector lambda; - - - double step_length,phi_zero,phi_alpha,phip,phip_zero; - double old_step,old_phi_zero,old_phip; - double L2_error; - double H1_error; - double Linfty_error; - double dual_L1_error; - double dual_L_infty_error; - FE_Q fe; - double p; - double line_search_tolerence; // c_1 in Nocedal & Wright - unsigned int dir_id; - std::string elements; - std::string Method; - -}; - -/*******************************************************************************************/ -// Boundary condition - - template - class BoundaryValues : public Function - { - public: - BoundaryValues () : Function() {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - }; - - - template - double BoundaryValues::value (const Point &pto, - const unsigned int /*component*/) const - { // could be anything else (theory works provided |Dg|_infty < 1/2) - return 0.0; - - /* A challenging BC leading to overdetermined problems - it is regulated by the parameter 00.5) & (theta<= pii-0.5)) - return eta*(theta-0.25); - else if ((theta>pii-0.5)&(theta<= pii+0.5)) - return eta*(pii-0.75-(theta-(pii-0.5))*(theta-(pii+0.5))); - else if ((theta>pii+0.5)&(theta<= 2*pii-0.5)) - return eta*((2*pii-theta)-0.25); - else - return eta*((theta-2*pii)*(theta-2*pii) );*/ - } - - - -/******************************************************************************/ -// Right-Hand Side -template -class RightHandSide : public Function -{ - public: - RightHandSide () : Function() {} - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - -template -double RightHandSide::value (const Point &p, -const unsigned int /*component*/) const -{ // set to constant = 4, for which explicit solution to compare exists - // could be anything - double return_value = 4.0; - return return_value; -} - - - -/*******************************************************************/ -// The ElastoplasticTorsion class implementation - - // Constructor of the class - template -ElastoplasticTorsion::ElastoplasticTorsion (ParameterHandler ¶m): - prm(param), - dof_handler (triangulation), - L2_error(1.0), - H1_error(1.0), - Linfty_error(1.0), - dual_L1_error(1.0), - dual_L_infty_error(1.0), - fe(2) - { - prm.enter_subsection ("Global Parameters"); - p=prm.get_double("p"); - prm.leave_subsection (); - prm.enter_subsection ("Algorithm Parameters"); - line_search_tolerence=prm.get_double("line_search_tolerence"); - dir_id=prm.get_integer("Descent_direction"); - prm.leave_subsection (); - if (fe.degree==1) - elements="P1"; - else elements="P2"; - - if (dir_id==0) - Method="Precond"; - else - Method="Newton"; -} - - - -template -ElastoplasticTorsion::~ElastoplasticTorsion () - { - dof_handler.clear (); - } - -/*****************************************************************************************/ -// print iteration message - -template -void ElastoplasticTorsion::print_it_message (const int counter, bool ks) -{ - if(ks){ - process_solution (counter); - std::cout << "iteration="<< counter+1 << " J(u_h)= "<< phi_zero << ", H1 error: " - << H1_error <<", W0-1,infty error: "<< Linfty_error<< " J'(u_h)(w)= "<< phip - << ", |J'(u_h)|= "<< system_rhs.l2_norm()< -void ElastoplasticTorsion::process_solution (const unsigned int it) -{ - Vector difference_per_cell (triangulation.n_active_cells()); - - // compute L2 error (save to difference_per_cell) - VectorTools::integrate_difference (dof_handler,present_solution, - Solution(),difference_per_cell,QGauss(3),VectorTools::L2_norm); - L2_error = difference_per_cell.l2_norm(); - - // compute H1 error (save to difference_per_cell) - VectorTools::integrate_difference (dof_handler,present_solution,Solution(), - difference_per_cell,QGauss(3),VectorTools::H1_seminorm); - H1_error = difference_per_cell.l2_norm(); - - // compute W1infty error (save to difference_per_cell) - const QTrapez<1> q_trapez; - const QIterated q_iterated (q_trapez, 5); - VectorTools::integrate_difference (dof_handler,present_solution,Solution(), - difference_per_cell,q_iterated,VectorTools::W1infty_seminorm); - Linfty_error = difference_per_cell.linfty_norm(); - - - convergence_table.add_value("cycle", it); - convergence_table.add_value("p", p); - convergence_table.add_value("L2", L2_error); - convergence_table.add_value("H1", H1_error); - convergence_table.add_value("Linfty", Linfty_error); - convergence_table.add_value("function value", phi_alpha); - convergence_table.add_value("derivative", phip); -} - - -/***************************************/ -// fill-in entry for the multiplier -template -void ElastoplasticTorsion::process_multiplier (const unsigned int cycle, const int iter,double time) -{ - const unsigned int n_active_cells=triangulation.n_active_cells(); - const unsigned int n_dofs=dof_handler.n_dofs(); - dual_L1_error=dual_error(); - dual_L_infty_error=dual_infty_error(); - - - dual_convergence_table.add_value("cycle", cycle); - dual_convergence_table.add_value("p", p); - dual_convergence_table.add_value("iteration_number", iter); - dual_convergence_table.add_value("cpu_time", time); - dual_convergence_table.add_value("cells", n_active_cells); - dual_convergence_table.add_value("dofs", n_dofs); - dual_convergence_table.add_value("L2", L2_error); - dual_convergence_table.add_value("H1", H1_error); - dual_convergence_table.add_value("Linfty", Linfty_error); - dual_convergence_table.add_value("dual_L1", dual_L1_error); - dual_convergence_table.add_value("dual_Linfty", dual_L_infty_error); - -} - - - - -/****************************************************************************************/ -// ElastoplasticTorsion::setup_system -// unchanged from step-15 - - template - void ElastoplasticTorsion::setup_system (const bool initial_step) - { - if (initial_step) - { - dof_handler.distribute_dofs (fe); - present_solution.reinit (dof_handler.n_dofs()); - grad_norm.reinit (dof_handler.n_dofs()); - lambda.reinit (dof_handler.n_dofs()); - - hanging_node_constraints.clear (); - DoFTools::make_hanging_node_constraints (dof_handler, - hanging_node_constraints); - hanging_node_constraints.close (); - } - - - // The remaining parts of the function - - 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); - } - -/***************************************************************************************/ - /* the coeffcients W, W' and G defining the problem. - - Min_u \int W(|Du|^2) dx - - They must be consistent as G(s)=W'(s)+2s W''(s) for any s>0. - recall that they receive the SQUARED gradient. */ - - template - double ElastoplasticTorsion::W (double Du2) const - { - return Du2; - } - - template - double ElastoplasticTorsion::Wp (double Du2) const - { - return 1.0; - } - - template - double ElastoplasticTorsion::G (double Du2) const - { - return 1.0; - } -/***************************************************************************************/ - - template - void ElastoplasticTorsion::assemble_system () - { - - - const QGauss quadrature_formula(3); - const RightHandSide right_hand_side; - system_matrix = 0; - system_rhs = 0; - - FEValues fe_values (fe, quadrature_formula, - update_gradients | - 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(); - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); - - std::vector > old_solution_gradients(n_q_points); - 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_matrix = 0; - cell_rhs = 0; - - fe_values.reinit (cell); - fe_values.get_function_gradients(present_solution, - old_solution_gradients); - - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) - { - long double coeff=0.0; - long double a=old_solution_gradients[q_point] * old_solution_gradients[q_point]; - long double exponent=(p-2.0)/2*std::log(a); - coeff= std::exp( exponent); - for (unsigned int i=0; iget_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); - } - - - - -/********************************** Refine Mesh ****************************************/ -// unchanged from step-15 - - template - void ElastoplasticTorsion::refine_mesh () - { - Vector estimated_error_per_cell (triangulation.n_active_cells()); - KellyErrorEstimator::estimate (dof_handler, - QGauss(3), - typename FunctionMap::type(), - present_solution, - estimated_error_per_cell); - - prm.enter_subsection ("Mesh & Refinement Parameters"); - const double top_fraction=prm.get_double("top_fraction_of_cells"); - const double bottom_fraction=prm.get_double("bottom_fraction_of_cells"); - prm.leave_subsection (); - GridRefinement::refine_and_coarsen_fixed_number (triangulation, - estimated_error_per_cell, - top_fraction, bottom_fraction); - - triangulation.prepare_coarsening_and_refinement (); - SolutionTransfer solution_transfer(dof_handler); - solution_transfer.prepare_for_coarsening_and_refinement(present_solution); - triangulation.execute_coarsening_and_refinement(); - dof_handler.distribute_dofs(fe); - Vector tmp(dof_handler.n_dofs()); - solution_transfer.interpolate(present_solution, tmp); - present_solution = tmp; - set_boundary_values (); - hanging_node_constraints.clear(); - - DoFTools::make_hanging_node_constraints(dof_handler, - hanging_node_constraints); - hanging_node_constraints.close(); - hanging_node_constraints.distribute (present_solution); - setup_system (false); - } - - -/*******************************************************************************************/ -// Dump the norm of the gradient and the lagrange multiplier in vtu format for visualization - template - void ElastoplasticTorsion::output_results (unsigned int counter) const - { - // multiplier object contains both |Du| and lambda. - ComputeMultiplier multiplier(p); - DataOut data_out; - - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (present_solution, "solution"); - data_out.add_data_vector (present_solution, multiplier); - data_out.build_patches (); - std::ostringstream p_str; - p_str << p<<"-cycle-"< - void ElastoplasticTorsion::set_boundary_values () - { - std::map boundary_values; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - BoundaryValues(), - boundary_values); - for (std::map::const_iterator - bp = boundary_values.begin(); - bp != boundary_values.end(); ++bp) - present_solution(bp->first) = bp->second; - } - - -/****************************************************************************************/ -// COMPUTE \phi(\alpha)=J_p(u_h+\alpha w) - template - double ElastoplasticTorsion::phi (const double alpha) const - { - double obj = 0.0; - const RightHandSide right_hand_side; - Vector evaluation_point (dof_handler.n_dofs()); - evaluation_point = present_solution; // copy of u_h - evaluation_point.add (alpha, newton_update); // u_{n+1}=u_n+alpha w_n - - const QGauss quadrature_formula(3); - FEValues fe_values (fe, quadrature_formula, - update_gradients | - 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(); - - Vector cell_residual (dofs_per_cell); - std::vector > gradients(n_q_points); - std::vector values(n_q_points); - - - 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_residual = 0; - fe_values.reinit (cell); - fe_values.get_function_gradients (evaluation_point, gradients); - fe_values.get_function_values (evaluation_point, values); - - - for (unsigned int q_point=0; q_point - double ElastoplasticTorsion::dual_error () const - { - double obj = 0.0; - - const QGauss quadrature_formula(3); - 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(); - - Vector cell_residual (dofs_per_cell); - std::vector > gradients(n_q_points); - - 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_residual = 0; - fe_values.reinit (cell); - fe_values.get_function_gradients (present_solution, gradients); - - for (unsigned int q_point=0; q_point0.5) - exact= 2*r-1; - - obj+=( std::abs(coeff-exact) ) * fe_values.JxW(q_point); - } - - } - - return obj; - } - -/*******************************************************************************************/ -// Compute L^infinity error norm of Lagrange Multiplier -// with respect to exact solution (cf. Alvarez & Flores, 2015) - - template - double ElastoplasticTorsion::dual_infty_error () const - { - double obj = 0.0; - const QTrapez<1> q_trapez; - const QIterated quadrature_formula (q_trapez, 10); - - FEValues fe_values (fe, quadrature_formula, - update_gradients | - update_quadrature_points ); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - Vector cell_residual (dofs_per_cell); - std::vector > gradients(n_q_points); - - 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_residual = 0; - fe_values.reinit (cell); - fe_values.get_function_gradients (present_solution, gradients); - - for (unsigned int q_point=0; q_point0.5) - exact= 2*r-1.0; - // compute |Du|^(p-2) as exp(p-2/2*log(Du^2)) - long double exponent=(p-2.0)/2*std::log(sqdGrad); - long double coeff=std::exp(exponent); - - if(std::abs(coeff-exact)>obj ) - obj=std::abs(coeff-exact); - } - - } - - return obj; - } - -/*****************************************************************************************/ -// check whether putative step-length satisfies sufficient decrease conditions - template -bool ElastoplasticTorsion::checkWolfe(double & alpha, double & phi_alpha) const -{ -if (phi_alpha< phi_zero+line_search_tolerence*phip*alpha ) - return true; - else - return false; -} - - -/*****************************************************************************************/ -// Find a step-length satisfying sufficient decrease condition by line-search -// uses quadratic interpolation - - template -bool ElastoplasticTorsion::determine_step_length(const int inner_it) -{ - unsigned int it=0; - bool done; - double alpha,nalpha; - prm.enter_subsection ("Algorithm Parameters"); - const unsigned int max_LS_it=prm.get_integer("max_LS_it"); - double init_SL=prm.get_double("init_step_length"); - prm.leave_subsection (); - if (inner_it==0) - alpha=init_SL; - else - { - alpha=std::min(1.45*old_step*old_phip/phip,1.0); - } - phi_alpha=phi(alpha); - std::cerr << "Step length=" << alpha << ", Value= " << phi_alpha; - // check if step-size satisfies sufficient decrease condition - done=checkWolfe(alpha,phi_alpha); - if (done) - std::cerr << " accepted" << std::endl; - else - std::cerr << " rejected" ; - - while ((!done) & (it1e3*std::abs(phi_zero) ) - nalpha=alpha/10; - alpha=nalpha; - phi_alpha=phi(alpha); - done=checkWolfe(alpha,phi_alpha); - if (done) - std::cerr << ", finished with steplength= "<< alpha<< ", fcn value= "<< phi_alpha< boundary; - triangulation.set_boundary (0, boundary); - } - else if (domain_id==1){ - // For the unit square - GridGenerator::hyper_cube (triangulation, 0, 1);} - else if (domain_id==2){ - /* For Glowinski's domain - ___ ___ __ 1 - | |__| | __ .8 - | | - | | - |__________| __ 0 - - | | | | - 0 .4 .6 1 - - */ - Triangulation tria1; - Triangulation tria2; - Triangulation tria3; - Triangulation tria4; - Triangulation tria5; - Triangulation tria6; - GridGenerator::hyper_rectangle(tria1, Point<2>(0.0,0.0), Point<2>(0.4,0.8)); - GridGenerator::hyper_rectangle(tria2, Point<2>(0.0,0.8), Point<2>(0.4,1.0)); - GridGenerator::hyper_rectangle(tria3, Point<2>(0.4,0.0), Point<2>(0.6,0.8)); - GridGenerator::hyper_rectangle(tria4, Point<2>(0.6,0.0), Point<2>(1.0,0.8)); - GridGenerator::hyper_rectangle(tria5, Point<2>(0.6,0.8), Point<2>(1.0,1.0)); - GridGenerator::merge_triangulations (tria1, tria2, tria6); - GridGenerator::merge_triangulations (tria6, tria3, tria6); - GridGenerator::merge_triangulations (tria6, tria4, tria6); - GridGenerator::merge_triangulations (tria6, tria5, triangulation); - } - // perform initial refinements - triangulation.refine_global(init_ref); -} - -/**************************************************************************************************/ - // ElastoplasticTorsion::solve(inner_it) - // Performs one inner iteration - - template - bool ElastoplasticTorsion::solve (const int inner_it) - { - prm.enter_subsection ("Algorithm Parameters"); - const unsigned int max_CG_it=prm.get_integer("Max_CG_it"); - const double CG_tol=prm.get_double("CG_tol"); - prm.leave_subsection (); - - SolverControl solver_control (max_CG_it,CG_tol); - SolverCG<> solver (solver_control); - - PreconditionSSOR<> preconditioner; - preconditioner.initialize(system_matrix,0.25); - - solver.solve (system_matrix, newton_update, system_rhs, - preconditioner); - hanging_node_constraints.distribute (newton_update); - /****** save current quantities for line-search **** */ - // Recall that phi(alpha)=J(u+alpha w) - old_step=step_length; - old_phi_zero=phi_zero; - phi_zero=phi(0); // phi(0)=J(u) - old_phip=phip; - phip=-1.0*(newton_update*system_rhs); //phi'(0)=J'(u) *w, rhs=-J'(u). - if (inner_it==0) - phip_zero=phip; - - if (phip>0){ // this should not happen, step back - std::cout << "Not a descent direction!" < -void ElastoplasticTorsion::run () -{ - - // get parameters - prm.enter_subsection ("Mesh & Refinement Parameters"); - const int adapt_ref=prm.get_integer("No of adaptive refinements"); - prm.leave_subsection (); - prm.enter_subsection ("Algorithm Parameters"); - const int max_inner=prm.get_integer("Max_inner"); - const double eps=prm.get_double("eps"); - const double hi_eps=prm.get_double("hi_eps"); - const int hi_th=prm.get_integer("hi_th"); - const double init_p=prm.get_double("init_p"); - const double delta_p=prm.get_double("delta_p"); - prm.leave_subsection (); - prm.enter_subsection ("Global Parameters"); - bool known_solution=prm.get_bool("known_solution"); - double actual_p=prm.get_double("p"); - prm.leave_subsection (); - /************************/ - - // init Timer - Timer timer; - double ptime=0.0; - timer.start (); - - // initalize mesh for the selected domain - init_mesh(); - - // setup FE space - setup_system (true); - set_boundary_values (); - - // init counters - int global_it=0; // Total inner iterations (counting both loops) - int cycle=0; // Total outer iterations (counting both loops) - int refinement = 0; // Refinements performed (adaptive) = outer iterations 2nd loop - - - // prepare to start first loop - p=init_p; - bool well_solved=true; - - /***************************** First loop ***********************************/ - /****************** Prepare initial condition using increasing p *************************/ - while(p=1)) | - !well_solved - ) - break; - } - ptime=timer(); - if (well_solved) - output_results (cycle); - - if(known_solution){ - process_multiplier(cycle,global_it,ptime); - //dual_convergence_table.write_tex(dual_error_table_file); - } - refine_mesh(); - cycle++; - p+=delta_p; - } - /*************************** first loop finished ********************/ - - - // prepare for second loop - p=actual_p; - well_solved=true; - - - /***************************** Second loop *********************************/ - /**************************** Solve problem for target p *********************************/ - - std::cout << "============ Solving problem with p=" <

ElastoplasticTorsionProblem(prm); - ElastoplasticTorsionProblem .run (); - } - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - return 0; -} diff --git a/Elastoplastic Torsion/doc/author b/Elastoplastic Torsion/doc/author deleted file mode 100644 index 8ff862d..0000000 --- a/Elastoplastic Torsion/doc/author +++ /dev/null @@ -1 +0,0 @@ -Salvador Flores diff --git a/Elastoplastic Torsion/doc/builds-on b/Elastoplastic Torsion/doc/builds-on deleted file mode 100644 index f54144f..0000000 --- a/Elastoplastic Torsion/doc/builds-on +++ /dev/null @@ -1 +0,0 @@ -step-15, step-29.