This program was contributed by Salvador Flores <sflores@dim.uchile.cl>.
It comes without any warranty or support by its authors or the authors of deal.II.
This program is part of the deal.II code gallery and consists of the following files (click to inspect):
Annotated version of ElastoplasticTorsion.cc
Include files
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/grid/
tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <typeinfo>
#include <fstream>
#include <iostream>
#include <deal.II/numerics/solution_transfer.h>
const ::Triangulation< dim, spacedim > & tria
Open a namespace for this program and import everything from the dealii namespace into it.
{
public:
void read_parameters(const std::string);
private:
void declare_parameters();
};
Constructor
prm(paramhandler)
{}
void ParameterReader::declare_parameters()
{
prm.enter_subsection ("Global Parameters");
{
"Penalization parameter");
"Whether the exact solution is known");
}
prm.leave_subsection ();
prm.enter_subsection ("Mesh & Refinement Parameters");
{
"Number identifying the domain in which we solve the problem");
"Number of global mesh refinement steps applied to initial coarse grid");
"Number of global adaptive mesh refinements");
"refinement threshold");
"coarsening threshold");
}
prm.leave_subsection ();
prm.enter_subsection ("Algorithm Parameters");
{
"0: Preconditioned descent, 1: Newton Method");
"Initial p");
"increase of p");
"Maximum Number of CG iterations");
"Tolerance for CG iterations");
"Maximum Number of LS iterations");
"line search tolerance constant (c1 in Nocedal-Wright)");
"initial step length in line-search");
"Maximum Number of inner iterations");
"Threshold on norm of the derivative to declare optimality achieved");
"Threshold on norm of the derivative to declare optimality achieved in highly refined mesh");
"Number of adaptive refinement before change convergence threshold");
}
prm.leave_subsection ();
}
void ParameterReader::read_parameters (const std::string parameter_file)
{
declare_parameters();
prm.parse_input (parameter_file);
}
The solution of the elastoplastic torsion problem on the unit disk with rhs=4.
template <int dim>
{
public:
virtual double value (
const Point<dim> &pto,
const unsigned int component = 0)
const override;
};
template <int dim>
double Solution<dim>::value (
const Point<dim> &pto,
const unsigned int)
const
{
double r=
sqrt(pto.square());
if (r<0.5)
else
return 1.0-r;
}
template <int dim>
{
double r=
sqrt(pto.square());
if (r<0.5)
return -2.0*pto;
else
return -1.0*pto/r;
}
virtual Tensor< 1, dim, RangeNumberType > gradient(const Point< dim > &p, const unsigned int component=0) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
template <int dim>
{
private:
double p;
public:
ComputeMultiplier (double pe);
virtual
void compute_derived_quantities_scalar (
const std::vector< double > &,
) const;
virtual std::vector<std::string>
get_names ()
const override;
virtual
std::vector<DataComponentInterpretation::DataComponentInterpretation>
};
template <int dim>
ComputeMultiplier<dim>::ComputeMultiplier (double pe): p(pe)
{}
void ComputeMultiplier<dim>::compute_derived_quantities_scalar(
const
std::vector< double > & ,
std::vector<
Vector< double > > &computed_quantities ) const
{
const unsigned
int n_quadrature_points = duh.size();
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
long double sqrGrad=duh[q]* duh[q];
long double exponent=(p-2.0)/2*
std::log(sqrGrad);
computed_quantities[q](0) =
std::sqrt(sqrGrad);
computed_quantities[q](1)=
std::exp(exponent);
}
}
template <int dim>
std::vector<std::string>
ComputeMultiplier<dim>::get_names() const
{
std::vector<std::string> solution_names;
solution_names.push_back ("Gradient norm");
solution_names.push_back ("Lagrange multiplier");
return solution_names;
}
template <int dim>
ComputeMultiplier<dim>::get_needed_update_flags () const
{
}
template <int dim>
std::vector<DataComponentInterpretation::DataComponentInterpretation>
ComputeMultiplier<dim>:: get_data_component_interpretation () const
{
std::vector<DataComponentInterpretation::DataComponentInterpretation>
interpretation;
virtual UpdateFlags get_needed_update_flags() const =0
virtual std::vector< std::string > get_names() const =0
virtual std::vector< DataComponentInterpretation::DataComponentInterpretation > get_data_component_interpretation() const
@ update_gradients
Shape function gradients.
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
norm of the gradient
Lagrange multiplier
template <int dim>
class ElastoplasticTorsion
{
public:
~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;
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;
double p;
double line_search_tolerence;
unsigned int dir_id;
std::string elements;
std::string Method;
};
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Boundary condition
template <int dim>
class BoundaryValues :
public Function<dim>
{
public:
const unsigned int component = 0) const override;
};
template <int dim>
double BoundaryValues<dim>::value (
const Point<dim> &,
const unsigned int ) const
{
could be anything else (theory works provided |Dg|_infty < 1/2)
Right-Hand Side
template <int dim>
class RightHandSide :
public Function<dim>
{
public:
const unsigned int component = 0) const override;
};
template <int dim>
double RightHandSide<dim>::value (
const Point<dim> &,
const unsigned int ) 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 <int dim>
prm(param),
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 <int dim>
ElastoplasticTorsion<dim>::~ElastoplasticTorsion ()
{
dof_handler.clear ();
}
print iteration message
template <int dim>
void ElastoplasticTorsion<dim>::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()<<std::endl;
}
else
{
std::cout << "iteration= " << counter+1 << " J(u_h)= "
<< phi_alpha << " J'(u_h)= "<< phip<<std::endl;
}
}
Convergence Tables
formating
template <int dim>
void ElastoplasticTorsion<dim>::format_convergence_tables()
{
convergence_table.set_precision("L2", 3);
convergence_table.set_precision("H1", 3);
convergence_table.set_precision("Linfty", 3);
convergence_table.set_precision("function value", 3);
convergence_table.set_precision("derivative", 3);
dual_convergence_table.set_precision("dual_L1", 3);
dual_convergence_table.set_precision("dual_Linfty", 3);
dual_convergence_table.set_precision("L2", 3);
dual_convergence_table.set_precision("H1", 3);
dual_convergence_table.set_precision("Linfty", 3);
convergence_table.set_scientific("L2", true);
convergence_table.set_scientific("H1", true);
convergence_table.set_scientific("Linfty", true);
convergence_table.set_scientific("function value", true);
convergence_table.set_scientific("derivative", true);
dual_convergence_table.set_scientific("dual_L1", true);
dual_convergence_table.set_scientific("dual_Linfty", true);
dual_convergence_table.set_scientific("L2", true);
dual_convergence_table.set_scientific("H1", true);
dual_convergence_table.set_scientific("Linfty", true);
}
fill-in entry for the solution
template <int dim>
void ElastoplasticTorsion<dim>::process_solution (const unsigned int it)
{
compute L2 error (save to difference_per_cell)
L2_error = difference_per_cell.l2_norm();
compute H1 error (save to difference_per_cell)
H1_error = difference_per_cell.l2_norm();
compute W1infty error (save to difference_per_cell)
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 <int dim>
void ElastoplasticTorsion<dim>::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 <int dim>
void ElastoplasticTorsion<dim>::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 ();
hanging_node_constraints);
hanging_node_constraints.close ();
}
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
The remaining parts of the function
newton_update.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
hanging_node_constraints.condense (c_sparsity);
sparsity_pattern.copy_from(c_sparsity);
system_matrix.reinit (sparsity_pattern);
}
template <int dim>
double ElastoplasticTorsion<dim>::W (double Du2) const
{
return Du2;
}
template <int dim>
double ElastoplasticTorsion<dim>::Wp (double ) const
{
return 1.0;
}
template <int dim>
double ElastoplasticTorsion<dim>::G (double ) const
{
return 1.0;
}
template <int dim>
void ElastoplasticTorsion<dim>::assemble_system ()
{
const RightHandSide<dim> right_hand_side;
system_matrix = 0;
system_rhs = 0;
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
std::vector<Tensor<1, dim> > old_solution_gradients(n_q_points);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
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);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
if (dir_id==1)
{
cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * fe_values.shape_grad(j, q_point)
* (G(a)+(p-1.0)*coeff) * fe_values.JxW(q_point);
}
else
{
cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * fe_values.shape_grad(j, q_point)
* (Wp(a)+coeff)
* fe_values.JxW(q_point);
}
}
cell_rhs(i) -= ( fe_values.shape_grad(i, q_point)
* old_solution_gradients[q_point]
* (Wp(a)+coeff)
-right_hand_side.value(fe_values.quadrature_point(q_point))
*fe_values.shape_value(i, q_point)
)
* fe_values.JxW(q_point);
}
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
hanging_node_constraints.condense (system_matrix);
hanging_node_constraints.condense (system_rhs);
std::map<types::global_dof_index,double> boundary_values;
0,
boundary_values);
system_matrix,
newton_update,
system_rhs);
}
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
unchanged from step-15
template <int dim>
void ElastoplasticTorsion<dim>::refine_mesh ()
{
using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
FunctionMap(),
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 ();
estimated_error_per_cell,
top_fraction, bottom_fraction);
solution_transfer.prepare_for_coarsening_and_refinement(present_solution);
dof_handler.distribute_dofs(fe);
solution_transfer.interpolate(present_solution, tmp);
present_solution = tmp;
set_boundary_values ();
hanging_node_constraints.clear();
hanging_node_constraints);
hanging_node_constraints.close();
hanging_node_constraints.distribute (present_solution);
setup_system (false);
}
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, typename InputVector::value_type > * > &neumann_bc, const InputVector &solution, Vector< float > &error, const ComponentMask &component_mask=ComponentMask(), const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
Dump the norm of the gradient and the lagrange multiplier in vtu format for visualization
template <int dim>
void ElastoplasticTorsion<dim>::output_results (unsigned int counter) const
{
multiplier object contains both |Du| and lambda.
ComputeMultiplier<dim> multiplier(p);
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-"<<counter;
std::string str = p_str.str();
const std::string filename = "solution-" + str+".vtu";
std::ofstream output (filename.c_str());
data_out.write_vtu (output);
}
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
unchanged from step-15
template <int dim>
void ElastoplasticTorsion<dim>::set_boundary_values ()
{
std::map<types::global_dof_index, double> boundary_values;
0,
BoundaryValues<dim>(),
boundary_values);
for (std::map<types::global_dof_index, double>::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 <int dim>
double ElastoplasticTorsion<dim>::phi (const double alpha) const
{
double obj = 0.0;
const RightHandSide<dim> right_hand_side;
evaluation_point = present_solution;
evaluation_point.add (alpha, newton_update);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
std::vector<Tensor<1, dim> > gradients(n_q_points);
std::vector<double> values(n_q_points);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
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<n_q_points; ++q_point)
{
double Du2=gradients[q_point] * gradients[q_point];
double penalty;
if (Du2<1.0e-10)
penalty=0.0;
else
obj+= 1/2 W(|Du|^2)+1/p |Du|^p -fu (see (1))
obj+=(
(0.5*W(Du2)+penalty/p)- right_hand_side.value(fe_values.quadrature_point(q_point))*values[q_point]
) * fe_values.JxW(q_point);
}
}
return obj;
}
Compute L^1 error norm of Lagrange Multiplier with respect to exact solution (cf. Alvarez & Flores, 2015)
template <int dim>
double ElastoplasticTorsion<dim>::dual_error () const
{
double obj = 0.0;
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
std::vector<Tensor<1, dim> > gradients(n_q_points);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
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_point<n_q_points; ++q_point)
{
double coeff=gradients[q_point] * gradients[q_point] ;
if (coeff<1.0e-15)
coeff=0.0;
else
double r=
std::sqrt(fe_values.quadrature_point(q_point).square());
double exact=0;
if (r>0.5)
exact= 2*r-1;
obj+=(
std::abs(coeff-exact) ) * fe_values.JxW(q_point);
}
}
return obj;
}
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Compute L^infinity error norm of Lagrange Multiplier with respect to exact solution (cf. Alvarez & Flores, 2015)
template <int dim>
double ElastoplasticTorsion<dim>::dual_infty_error () const
{
double obj = 0.0;
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
std::vector<Tensor<1, dim> > gradients(n_q_points);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
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_point<n_q_points; ++q_point)
{
long double sqdGrad=gradients[q_point] * gradients[q_point] ;
double r=
std::sqrt(fe_values.quadrature_point(q_point).square());
double exact=0;
if (r>0.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);
}
}
return obj;
}
check whether putative step-length satisfies sufficient decrease conditions
template <int dim>
bool ElastoplasticTorsion<dim>::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 <int dim>
bool ElastoplasticTorsion<dim>::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;
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
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) & (it<max_LS_it))
{
new try obtained by quadratic interpolation
nalpha=-(phip*alpha*alpha)/(2*(phi_alpha-phi_zero-phip*alpha));
if (nalpha<1e-3*alpha ||
std::abs(nalpha-alpha)/alpha<1e-8)
nalpha=alpha/2;
else if ( phi_alpha-phi_zero>1e3*
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<<std::endl;
it=it+1;
}
if (!done)
{
std::cerr << ", max. no. of iterations reached wiht steplength= "<< alpha
<< ", fcn value= "<< phi_alpha<<std::endl;
return false;
}
else
{
step_length=alpha;
return true;
}
}
ElastoplasticTorsion::init_mesh()
template <int dim>
void ElastoplasticTorsion<dim>::init_mesh ()
{
get parameters
prm.enter_subsection ("Mesh & Refinement Parameters");
const int domain_id=prm.get_integer("Code for the domain");
const int init_ref=prm.get_integer("No of initial refinements");
prm.leave_subsection ();
if (domain_id==0)
{
For the unit disk around the origin
}
else if (domain_id==1)
{
void hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
For the unit square
}
else if (domain_id==2)
{
}
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
perform initial refinements
ElastoplasticTorsion::solve(inner_it) Performs one inner iteration
template <int dim>
bool ElastoplasticTorsion<dim>::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 ();
solver.solve (system_matrix, newton_update, system_rhs,
preconditioner);
hanging_node_constraints.distribute (newton_update);
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
Recall that phi(alpha)=J(u+alpha w)
old_step=step_length;
old_phi_zero=phi_zero;
phi_zero=phi(0);
old_phip=phip;
phip=-1.0*(newton_update*system_rhs);
if (inner_it==0)
phip_zero=phip;
if (phip>0)
{
std::cout << "Not a descent direction!" <<std::endl;
present_solution.add (-1.0*step_length, newton_update);
step_length=step_length/2;
phip=old_phip;
return false;
}
else
{
if (determine_step_length(inner_it))
{
update u_{n+1}=u_n+alpha w_n
present_solution.add (step_length, newton_update);
return true;
}
else return false;
}
}
ElastoplasticTorsion::run
template <int dim>
void ElastoplasticTorsion<dim>::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
initalize mesh for the selected domain
setup FE space
setup_system (true);
set_boundary_values ();
init counters
int global_it=0;
int cycle=0;
int refinement = 0;
prepare to start first loop
p=init_p;
bool well_solved=true;
while (p<actual_p)
{
std::cout <<"--Preparing initial condition with p="<<p<<" iter.= " << global_it<< " .-- "<< std::endl;
timer.restart();
for (int inner_iteration=0; inner_iteration<max_inner; ++inner_iteration,++global_it)
{
assemble_system ();
well_solved=solve (inner_iteration);
print_it_message (global_it, known_solution);
if (
((system_rhs.l2_norm()/
std::sqrt(system_rhs.size()) <1e-4) & (cycle<1)) |
((system_rhs.l2_norm()/
std::sqrt(system_rhs.size()) <1e-5) & (cycle>=1)) |
!well_solved
)
break;
}
ptime=timer.cpu_time();
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;
}
prepare for second loop
p=actual_p;
well_solved=true;
std::cout << "============ Solving problem with p=" <<p << " ==================" << std::endl;
while ((cycle<adapt_ref) & well_solved)
{
timer.restart();
inner iteration
for (int inner_iteration=0; inner_iteration<max_inner; ++inner_iteration,++global_it)
{
assemble_system ();
well_solved=solve (inner_iteration);
print_it_message (global_it, known_solution);
if (
((system_rhs.l2_norm()/
std::sqrt(system_rhs.size()) < eps) & (refinement<hi_th)) |
(( system_rhs.l2_norm()/
std::sqrt (system_rhs.size()) <hi_eps) | (!well_solved))
)
break;
}
inner iterations finished
ptime=timer.cpu_time();
if (well_solved)
output_results (cycle);
compute and display error, if the explicit solution is known
if (known_solution)
{
process_multiplier(cycle,global_it,ptime);
std::cout << "finished with H1 error: " << H1_error << ", dual error (L1): "
<< dual_L1_error << "dual error (L infty): "<<dual_L_infty_error <<std::endl;
}
update counters
refine mesh
std::cout << "******** Refined mesh " << cycle << " ********" << std::endl;
refine_mesh();
}
write convergence tables to file
if (known_solution)
{
format_convergence_tables();
std::string error_filename = "error"+Method+elements+".tex";
std::ofstream error_table_file(error_filename.c_str());
std::string dual_error_filename = "dual_error"+Method+elements+".tex";
std::ofstream dual_error_table_file(dual_error_filename.c_str());
convergence_table.write_tex(error_table_file);
dual_convergence_table.write_tex(dual_error_table_file);
}
}
}
The main function
int main ()
{
try
{
using namespace nsp;
ParameterReader param(prm);
param.read_parameters("EPT.prm");
ElastoplasticTorsion<2> 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;
}
unsigned int depth_console(const unsigned int n)