From d58329080d6370eaed4401c0768f28430b7c8fbe Mon Sep 17 00:00:00 2001 From: Salvador Flores Date: Fri, 16 Oct 2015 16:02:26 -0300 Subject: [PATCH] Solve p-laplacian elliptic problems for large p --- Elastoplastic Torsion/CMakeLists.txt | 40 + Elastoplastic Torsion/ElastoplasticTorsion.cc | 1343 +++++++++++++++++ Elastoplastic Torsion/doc/author | 1 + Elastoplastic Torsion/doc/builds-on | 1 + 4 files changed, 1385 insertions(+) create mode 100644 Elastoplastic Torsion/CMakeLists.txt create mode 100644 Elastoplastic Torsion/ElastoplasticTorsion.cc create mode 100644 Elastoplastic Torsion/doc/author create mode 100644 Elastoplastic Torsion/doc/builds-on diff --git a/Elastoplastic Torsion/CMakeLists.txt b/Elastoplastic Torsion/CMakeLists.txt new file mode 100644 index 0000000..82fb113 --- /dev/null +++ b/Elastoplastic Torsion/CMakeLists.txt @@ -0,0 +1,40 @@ +## +# 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 new file mode 100644 index 0000000..1cec33a --- /dev/null +++ b/Elastoplastic Torsion/ElastoplasticTorsion.cc @@ -0,0 +1,1343 @@ +/* --------------------------------------------------------------------- + * + * 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 new file mode 100644 index 0000000..8ff862d --- /dev/null +++ b/Elastoplastic Torsion/doc/author @@ -0,0 +1 @@ +Salvador Flores diff --git a/Elastoplastic Torsion/doc/builds-on b/Elastoplastic Torsion/doc/builds-on new file mode 100644 index 0000000..f54144f --- /dev/null +++ b/Elastoplastic Torsion/doc/builds-on @@ -0,0 +1 @@ +step-15, step-29. -- 2.39.5