--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * 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
+ "Solving variational problems with uniform gradient bounds by p-Laplacian
+ approximation: Elastoplastic torsion implementation using the deal.II
+ library"
+ available at the documentation and at http://www.dim.uchile.cl/~sflores.
+
+*/
+
+// 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/compressed_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/constraint_matrix.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/tria_boundary_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>
+
+// 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 <int dim>
+ class Solution : public Function<dim>
+ {
+ public:
+ Solution () : Function<dim>() {}
+ virtual double value (const Point<dim> &pto, const unsigned int component = 0) const;
+ virtual Tensor<1,dim> gradient (const Point<dim> &pto, const unsigned int component = 0) const;
+ };
+
+ template <int dim>
+ double Solution<dim>::value (const Point<dim> &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 <int dim>
+ Tensor<1,dim> Solution<dim>::gradient (const Point<dim> &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 <int dim>
+ class ComputeMultiplier : public DataPostprocessor<dim>
+ {
+ 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<std::string> get_names () const;
+
+ virtual
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ get_data_component_interpretation () const;
+ virtual UpdateFlags get_needed_update_flags () const;
+ };
+
+
+ template <int dim>
+ ComputeMultiplier<dim>::ComputeMultiplier (double pe): p(pe)
+ {}
+
+
+ template <int dim>
+ void ComputeMultiplier<dim>::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<n_quadrature_points; ++q)
+ {
+ long double sqrGrad=duh[q]* duh[q]; //squared norm of the gradient
+ long double exponent=(p-2.0)/2*std::log(sqrGrad);
+ computed_quantities[q](0) = std::sqrt(sqrGrad); // norm of the gradient
+ computed_quantities[q](1)= std::exp(exponent); // multiplier
+ }
+ }
+
+
+
+
+
+ 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>
+ UpdateFlags
+ ComputeMultiplier<dim>::get_needed_update_flags () const
+ {
+ return update_gradients;
+ }
+
+
+
+ template <int dim>
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ ComputeMultiplier<dim>:: get_data_component_interpretation () const
+ {
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ interpretation;
+ // norm of the gradient
+ interpretation.push_back (DataComponentInterpretation::component_is_scalar);
+ // Lagrange multiplier
+ interpretation.push_back (DataComponentInterpretation::component_is_scalar);
+ return interpretation;
+ }
+
+
+
+
+
+// *************************************************************************************** //
+ template <int dim>
+ 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<dim> triangulation;
+ DoFHandler<dim> dof_handler;
+ ConstraintMatrix hanging_node_constraints;
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+ ConvergenceTable convergence_table;
+ ConvergenceTable dual_convergence_table;
+ Vector<double> present_solution;
+ Vector<double> newton_update;
+ Vector<double> system_rhs;
+ Vector<double> grad_norm;
+ Vector<double> 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<dim> 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 <int dim>
+ class BoundaryValues : public Function<dim>
+ {
+ public:
+ BoundaryValues () : Function<dim>() {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+ };
+
+
+ template <int dim>
+ double BoundaryValues<dim>::value (const Point<dim> &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 0<eta<1.
+ eta closer to 1 leads to more difficult problems.
+
+ double pii=numbers::PI;
+ double theta=std::atan2(p[1],p[0])+pii;
+ double eta=0.9;
+
+ if (theta <= 0.5)
+ return eta*(theta*theta);
+ else if ((theta >0.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 <int dim>
+ class RightHandSide : public Function<dim>
+ {
+ public:
+ RightHandSide () : Function<dim>() {}
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+ };
+
+ template <int dim>
+ double RightHandSide<dim>::value (const Point<dim> &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 <int dim>
+ ElastoplasticTorsion<dim>::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 <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)
+ {
+ Vector<float> difference_per_cell (triangulation.n_active_cells());
+
+ // compute L2 error (save to difference_per_cell)
+ VectorTools::integrate_difference (dof_handler,present_solution,
+ Solution<dim>(),difference_per_cell,QGauss<dim>(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<dim>(),
+ difference_per_cell,QGauss<dim>(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<dim> q_iterated (q_trapez, 5);
+ VectorTools::integrate_difference (dof_handler,present_solution,Solution<dim>(),
+ 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 <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 ();
+ 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 <int dim>
+ double ElastoplasticTorsion<dim>::W (double Du2) const
+ {
+ return Du2;
+ }
+
+ template <int dim>
+ double ElastoplasticTorsion<dim>::Wp (double Du2) const
+ {
+ return 1.0;
+ }
+
+ template <int dim>
+ double ElastoplasticTorsion<dim>::G (double Du2) const
+ {
+ return 1.0;
+ }
+ /***************************************************************************************/
+
+ template <int dim>
+ void ElastoplasticTorsion<dim>::assemble_system ()
+ {
+
+
+ const QGauss<dim> quadrature_formula(3);
+ const RightHandSide<dim> right_hand_side;
+ system_matrix = 0;
+ system_rhs = 0;
+
+ FEValues<dim> 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<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+
+ std::vector<Tensor<1, dim> > old_solution_gradients(n_q_points);
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+
+ typename DoFHandler<dim>::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; 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],
+ cell_matrix(i,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;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ZeroFunction<dim>(),
+ boundary_values);
+ MatrixTools::apply_boundary_values (boundary_values,
+ system_matrix,
+ newton_update,
+ system_rhs);
+ }
+
+
+
+
+ /********************************** Refine Mesh ****************************************/
+// unchanged from step-15
+
+ template <int dim>
+ void ElastoplasticTorsion<dim>::refine_mesh ()
+ {
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ QGauss<dim-1>(3),
+ typename FunctionMap<dim>::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<dim> solution_transfer(dof_handler);
+ solution_transfer.prepare_for_coarsening_and_refinement(present_solution);
+ triangulation.execute_coarsening_and_refinement();
+ dof_handler.distribute_dofs(fe);
+ Vector<double> 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 <int dim>
+ void ElastoplasticTorsion<dim>::output_results (unsigned int counter) const
+ {
+ // multiplier object contains both |Du| and lambda.
+ ComputeMultiplier<dim> multiplier(p);
+ DataOut<dim> 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-"<<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);
+ }
+
+ /********************************************************************************************/
+// unchanged from step-15
+ template <int dim>
+ void ElastoplasticTorsion<dim>::set_boundary_values ()
+ {
+ std::map<types::global_dof_index, double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 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;
+ Vector<double> 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<dim> quadrature_formula(3);
+ FEValues<dim> 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<double> cell_residual (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);
+
+ typename DoFHandler<dim>::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<n_q_points; ++q_point)
+ {
+ double Du2=gradients[q_point] * gradients[q_point]; // Du2=|Du|^2
+ double penalty;
+ if (Du2<1.0e-10)
+ penalty=0.0;
+ else
+ penalty=std::pow(Du2,p/2.0); // penalty=|Du|^p
+
+ // 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 QGauss<dim> quadrature_formula(3);
+ FEValues<dim> 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<double> cell_residual (dofs_per_cell);
+ std::vector<Tensor<1, dim> > gradients(n_q_points);
+
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::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_point<n_q_points; ++q_point)
+ {
+ double coeff=gradients[q_point] * gradients[q_point] ;
+ if (coeff<1.0e-15)
+ coeff=0.0;
+ else
+ coeff=std::pow(coeff,(p-2.0)/2.0); // |Du_p|^(p-2)
+
+ 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;
+ }
+
+ /*******************************************************************************************/
+// 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 QTrapez<1> q_trapez;
+ const QIterated<dim> quadrature_formula (q_trapez, 10);
+
+ FEValues<dim> 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<double> cell_residual (dofs_per_cell);
+ std::vector<Tensor<1, dim> > gradients(n_q_points);
+
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::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_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);
+ 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 <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;
+ // 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
+ GridGenerator::hyper_ball (triangulation);
+ static const HyperBallBoundary<dim> 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<dim> tria1;
+ Triangulation<dim> tria2;
+ Triangulation<dim> tria3;
+ Triangulation<dim> tria4;
+ Triangulation<dim> tria5;
+ Triangulation<dim> 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 <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 ();
+
+ 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!" <<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
+ 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<actual_p) // outer iteration, increasing 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();
+ 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=" <<p << " ==================" << std::endl;
+ /***** Outer iteration - refining mesh ******************/
+ 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();
+ 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
+ ++refinement;
+ ++cycle;
+ // refine mesh
+ std::cout << "******** Refined mesh " << cycle << " ********" << std::endl;
+ refine_mesh();
+ }// second loop
+
+ // 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);
+ }
+ }//run()
+
+}//namespace
+
+/**********************************************************************************************/
+// The main function
+int main ()
+{
+ try
+ {
+ using namespace dealii;
+ using namespace nsp;
+ deallog.depth_console (0);
+
+ ParameterHandler prm;
+ 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;
+}
projects / code-gallery.git / commitdiff