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):
Include files
Open a namespace for this program and import everything from the dealii namespace into it.
void declare_parameters();
Constructor
prm.enter_subsection (
"Global Parameters");
"Penalization parameter");
"Whether the exact solution is known");
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");
prm.enter_subsection (
"Algorithm Parameters");
"0: Preconditioned descent, 1: Newton Method");
"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");
{
The solution of the elastoplastic torsion problem on the unit disk with rhs=4.
virtual double value (
const Point<dim> &
pto,
const unsigned int component = 0)
const override;
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const std::vector< double > &,
virtual std::vector<std::string> get_names ()
const override;
std::vector<DataComponentInterpretation::DataComponentInterpretation>
get_data_component_interpretation ()
const override;
virtual UpdateFlags get_needed_update_flags ()
const override;
for (
unsigned int q=0;
q<n_quadrature_points; ++
q)
std::vector<DataComponentInterpretation::DataComponentInterpretation>
std::vector<DataComponentInterpretation::DataComponentInterpretation>
@ 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
double W (
double Du2)
const;
double Wp (
double Du2)
const;
double G (
double Du2)
const;
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Boundary condition
const unsigned int component = 0)
const override;
const unsigned int )
const
could be anything else (theory works provided |Dg|_infty < 1/2)
Right-Hand Side
const unsigned int component = 0)
const override;
const unsigned int )
const
set to constant = 4, for which explicit solution to compare exists could be anything
double return_value = 4.0;
The ElastoplasticTorsion class implementation
Constructor of the class
prm.enter_subsection (
"Global Parameters");
prm.enter_subsection (
"Algorithm Parameters");
dir_id=prm.get_integer(
"Descent_direction");
print iteration message
{
std::cout <<
"iteration="<< counter+1 <<
" J(u_h)= "<<
phi_zero <<
", H1 error: "
<<
", |J'(u_h)|= "<<
system_rhs.l2_norm()<<std::endl;
std::cout <<
"iteration= " << counter+1 <<
" J(u_h)= "
Convergence Tables
formatting
fill-in entry for the solution
compute L2 error (save to difference_per_cell)
compute H1 error (save to difference_per_cell)
compute W1infty error (save to difference_per_cell)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
fill-in entry for the multiplier
{
const unsigned int n_active_cells=
triangulation.n_active_cells();
const unsigned int n_dofs=dof_handler.n_dofs();
ElastoplasticTorsion::setup_system unchanged from step-15
{
dof_handler.distribute_dofs (fe);
lambda.reinit (dof_handler.n_dofs());
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
The remaining parts of the function
system_matrix.reinit (sparsity_pattern);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
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)
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)
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],
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={}, 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
using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
prm.enter_subsection (
"Mesh & Refinement Parameters");
const double top_fraction=prm.get_double(
"top_fraction_of_cells");
dof_handler.distribute_dofs(fe);
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, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, 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)
#define DEAL_II_VERSION_GTE(major, minor, subminor)
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
multiplier object contains both |Du| and lambda.
data_out.attach_dof_handler (dof_handler);
data_out.build_patches ();
std::ostringstream
p_str;
p_str << p<<
"-cycle-"<<counter;
std::ofstream output (
filename.c_str());
data_out.write_vtu (output);
unchanged from step-15
for (std::map<types::global_dof_index, double>::const_iterator
COMPUTE \(\phi(\alpha)=J_p(u_h+\alpha w)\)
const unsigned int dofs_per_cell = fe.dofs_per_cell;
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)
obj+= 1/2 W(|Du|^2)+1/p |Du|^p -fu (see (1))
Compute L^1 error norm of Lagrange Multiplier with respect to exact solution (cf. Alvarez & Flores, 2015)
const unsigned int dofs_per_cell = fe.dofs_per_cell;
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)
::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)
const unsigned int dofs_per_cell = fe.dofs_per_cell;
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)
compute |Du|^(p-2) as exp(p-2/2*log(Du^2))
check whether putative step-length satisfies sufficient decrease conditions
Find a step-length satisfying sufficient decrease condition by line-search uses quadratic interpolation
{
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");
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
check if step-size satisfies sufficient decrease condition
std::cerr <<
" accepted" << std::endl;
std::cerr <<
" rejected" ;
new try obtained by quadratic interpolation
std::cerr <<
", finished with steplength= "<<
alpha<<
", fcn value= "<<
phi_alpha<<std::endl;
std::cerr <<
", max. no. of iterations reached with steplength= "<<
alpha
ElastoplasticTorsion::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");
For the unit disk around the origin
void hyper_ball(Triangulation< dim, spacedim > &tria, const Point< spacedim > ¢er={}, const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
For the unit square
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
{
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");
preconditioner.initialize(system_matrix,0.25);
Recall that phi(alpha)=J(u+alpha w)
std::cout <<
"Not a descent direction!" <<std::endl;
update u_{n+1}=u_n+alpha w_n
ElastoplasticTorsion::run
get parameters
prm.enter_subsection (
"Mesh & Refinement Parameters");
const int adapt_ref=prm.get_integer(
"No of adaptive refinements");
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.enter_subsection (
"Global Parameters");
init Timer
initialize mesh for the selected domain
setup FE space
init counters
prepare to start first loop
std::cout <<
"--Preparing initial condition with p="<<p<<
" iter.= " <<
global_it<<
" .-- "<< std::endl;
dual_convergence_table.write_tex(dual_error_table_file);
prepare for second loop
std::cout <<
"============ Solving problem with p=" <<p <<
" ==================" << std::endl;
inner iteration
inner iterations finished
compute and display error, if the explicit solution is known
std::cout <<
"finished with H1 error: " <<
H1_error <<
", dual error (L1): "
update counters
refine mesh
std::cout <<
"******** Refined mesh " << cycle <<
" ********" << std::endl;
write convergence tables to file
The main function
param.read_parameters(
"EPT.prm");
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::cerr << std::endl << std::endl
<<
"----------------------------------------------------"
std::cerr <<
"Unknown exception!" << std::endl
<<
"Aborting!" << std::endl
<<
"----------------------------------------------------"
unsigned int depth_console(const unsigned int n)