Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
ElastoplasticTorsion.h
Go to the documentation of this file.
1 
48  *
49  *
50  * /*
51  * This piece of software solves the elliptic p-laplacian
52  * boundary-value problems:
53  *
54  * Min {∫ 1/2 W(|Du|²)+ 1/p |Du|^p -fu : u=g on ∂S } (1)
55  * u
56  *
57  * for large values of p, which approximates (see Alvarez & Flores 2015)
58  *
59  * Min {∫ 1/2 W(|Du|²) -fu : |Du|<1 a.s. on S, u=g on ∂S }
60  * u
61  *
62  * By default W(t)=t and S=unit disk.
63  *
64  * Large portions of this code are borrowed from the deal.ii tutorials
65  *
66  * @ref step_15 "step-15" @ref step_29 "step-29".
67  *
68  * For further details see the technical report
69  * "Solving variational problems with uniform gradient bounds by p-Laplacian
70  * approximation: Elastoplastic torsion implementation using the deal.II
71  * library"
72  * available at the documentation and at http://www.dim.uchile.cl/~sflores.
73  *
74  * */
75  *
76  * @endcode
77  *
78  * Include files
79  *
80 
81  *
82  *
83  * @code
84  * #include <deal.II/base/quadrature_lib.h>
85  * #include <deal.II/base/function.h>
86  * #include <deal.II/base/logstream.h>
87  * #include <deal.II/base/utilities.h>
88  * #include <deal.II/base/convergence_table.h>
89  * #include <deal.II/base/smartpointer.h>
90  * #include <deal.II/base/parameter_handler.h>
91  * #include <deal.II/base/timer.h>
92  *
93  * #include <deal.II/lac/vector.h>
94  * #include <deal.II/lac/full_matrix.h>
95  * #include <deal.II/lac/sparse_matrix.h>
96  * #include <deal.II/lac/solver_cg.h>
97  * #include <deal.II/lac/precondition.h>
98  * #include <deal.II/lac/affine_constraints.h>
99  * #include <deal.II/lac/dynamic_sparsity_pattern.h>
100  *
101  * #include <deal.II/grid/tria.h>
102  * #include <deal.II/grid/grid_generator.h>
103  * #include <deal.II/grid/tria_accessor.h>
104  * #include <deal.II/grid/tria_iterator.h>
105  * #include <deal.II/grid/manifold_lib.h>
106  * #include <deal.II/grid/grid_refinement.h>
107  * #include <deal.II/grid/grid_in.h>
108  *
109  * #include <deal.II/dofs/dof_handler.h>
110  * #include <deal.II/dofs/dof_accessor.h>
111  * #include <deal.II/dofs/dof_tools.h>
112  * #include <deal.II/dofs/dof_renumbering.h>
113  *
114  * #include <deal.II/fe/fe_values.h>
115  * #include <deal.II/fe/fe_q.h>
116  *
117  * #include <deal.II/numerics/vector_tools.h>
118  * #include <deal.II/numerics/matrix_tools.h>
119  * #include <deal.II/numerics/data_out.h>
120  * #include <deal.II/numerics/error_estimator.h>
121  * #include <deal.II/numerics/solution_transfer.h>
122  *
123  * #include <typeinfo>
124  * #include <fstream>
125  * #include <iostream>
126  *
127  * #include <deal.II/numerics/solution_transfer.h>
128  *
129  * @endcode
130  *
131  * Open a namespace for this program and import everything from the
132  * dealii namespace into it.
133  *
134  * @code
135  * namespace nsp
136  * {
137  * using namespace dealii;
138  *
139  * @endcode
140  *
141  * ********************************************************
142  *
143  * @code
144  * class ParameterReader : public Subscriptor
145  * {
146  * public:
147  * ParameterReader(ParameterHandler &);
148  * void read_parameters(const std::string);
149  * private:
150  * void declare_parameters();
151  * ParameterHandler &prm;
152  * };
153  * @endcode
154  *
155  * Constructor
156  *
157  * @code
158  * ParameterReader::ParameterReader(ParameterHandler &paramhandler):
159  * prm(paramhandler)
160  * {}
161  *
162  * void ParameterReader::declare_parameters()
163  * {
164  *
165  * prm.enter_subsection ("Global Parameters");
166  * {
167  * prm.declare_entry("p", "100",Patterns::Double(2.1),
168  * "Penalization parameter");
169  * prm.declare_entry("known_solution", "true",Patterns::Bool(),
170  * "Whether the exact solution is known");
171  * }
172  * prm.leave_subsection ();
173  *
174  * prm.enter_subsection ("Mesh & Refinement Parameters");
175  * {
176  * prm.declare_entry("Code for the domain", "0",Patterns::Integer(0,2),
177  * "Number identifying the domain in which we solve the problem");
178  * prm.declare_entry("No of initial refinements", "4",Patterns::Integer(0),
179  * "Number of global mesh refinement steps applied to initial coarse grid");
180  * prm.declare_entry("No of adaptive refinements", "8",Patterns::Integer(0),
181  * "Number of global adaptive mesh refinements");
182  * prm.declare_entry("top_fraction_of_cells", "0.25",Patterns::Double(0),
183  * "refinement threshold");
184  * prm.declare_entry("bottom_fraction_of_cells", "0.05",Patterns::Double(0),
185  * "coarsening threshold");
186  * }
187  * prm.leave_subsection ();
188  *
189  *
190  * prm.enter_subsection ("Algorithm Parameters");
191  * {
192  * prm.declare_entry("Descent_direction", "0",Patterns::Integer(0,1),
193  * "0: Preconditioned descent, 1: Newton Method");
194  * prm.declare_entry("init_p", "10",Patterns::Double(2),
195  * "Initial p");
196  * prm.declare_entry("delta_p", "50",Patterns::Double(0),
197  * "increase of p");
198  * prm.declare_entry("Max_CG_it", "1500",Patterns::Integer(1),
199  * "Maximum Number of CG iterations");
200  * prm.declare_entry("CG_tol", "1e-10",Patterns::Double(0),
201  * "Tolerance for CG iterations");
202  * prm.declare_entry("max_LS_it", "45",Patterns::Integer(1),
203  * "Maximum Number of LS iterations");
204  * prm.declare_entry("line_search_tolerence", "1e-6",Patterns::Double(0),
205  * "line search tolerance constant (c1 in Nocedal-Wright)");
206  * prm.declare_entry("init_step_length", "1e-2",Patterns::Double(0),
207  * "initial step length in line-search");
208  * prm.declare_entry("Max_inner", "800",Patterns::Integer(1),
209  * "Maximum Number of inner iterations");
210  * prm.declare_entry("eps", "1.0e-8",Patterns::Double(0),
211  * "Threshold on norm of the derivative to declare optimality achieved");
212  * prm.declare_entry("hi_eps", "1.0e-9",Patterns::Double(0),
213  * "Threshold on norm of the derivative to declare optimality achieved in highly refined mesh");
214  * prm.declare_entry("hi_th", "8",Patterns::Integer(0),
215  * "Number of adaptive refinement before change convergence threshold");
216  * }
217  * prm.leave_subsection ();
218  *
219  * }
220  * void ParameterReader::read_parameters (const std::string parameter_file)
221  * {
222  * declare_parameters();
223  * prm.parse_input (parameter_file);
224  * }
225  *
226  * @endcode
227  *
228  * ******************************************************************************************
229  * The solution of the elastoplastic torsion problem on the unit disk with rhs=4.
230  *
231 
232  *
233  *
234  * @code
235  * template <int dim>
236  * class Solution : public Function<dim>
237  * {
238  * public:
239  * Solution () : Function<dim>() {}
240  * virtual double value (const Point<dim> &pto, const unsigned int component = 0) const;
241  * virtual Tensor<1,dim> gradient (const Point<dim> &pto, const unsigned int component = 0) const;
242  * };
243  *
244  * template <int dim>
245  * double Solution<dim>::value (const Point<dim> &pto,const unsigned int) const
246  * {
247  * double r=sqrt(pto.square());
248  * if (r<0.5)
249  * return -1.0*std::pow(r,2.0)+0.75;
250  * else
251  * return 1.0-r;
252  * }
253  *
254  *
255  *
256  * template <int dim>
257  * Tensor<1,dim> Solution<dim>::gradient (const Point<dim> &pto,const unsigned int) const
258  * {
259  * double r=sqrt(pto.square());
260  * if (r<0.5)
261  * return -2.0*pto;
262  * else
263  * return -1.0*pto/r;
264  * }
265  *
266  *
267  *
268  *
269  * @endcode
270  *
271  * ******************************************************************************************
272  *
273  * @code
274  * /* Compute the Lagrange multiplier (as a derived quantity) */
275  *
276  *
277  * template <int dim>
278  * class ComputeMultiplier : public DataPostprocessor<dim>
279  * {
280  * private:
281  * double p;
282  * public:
283  * ComputeMultiplier (double pe);
284  *
285  * virtual
286  * void compute_derived_quantities_scalar (
287  * const std::vector< double > &,
288  * const std::vector< Tensor< 1, dim > > &,
289  * const std::vector< Tensor< 2, dim > > &,
290  * const std::vector< Point< dim > > &,
291  * const std::vector< Point< dim > > &,
292  * std::vector< Vector< double > > &
293  * ) const;
294  *
295  * virtual std::vector<std::string> get_names () const;
296  *
297  * virtual
298  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
300  * virtual UpdateFlags get_needed_update_flags () const;
301  * };
302  *
303  *
304  * template <int dim>
305  * ComputeMultiplier<dim>::ComputeMultiplier (double pe): p(pe)
306  * {}
307  *
308  *
309  * template <int dim>
310  * void ComputeMultiplier<dim>::compute_derived_quantities_scalar(
311  * const std::vector< double > & /*uh*/,
312  * const std::vector< Tensor< 1, dim > > &duh,
313  * const std::vector< Tensor< 2, dim > > & /*dduh*/,
314  * const std::vector< Point< dim > > & /* normals*/,
315  * const std::vector< Point< dim > > & /*evaluation_points*/,
316  * std::vector< Vector< double > > &computed_quantities ) const
317  * {
318  * const unsigned int n_quadrature_points = duh.size();
319  *
320  * for (unsigned int q=0; q<n_quadrature_points; ++q)
321  * {
322  * long double sqrGrad=duh[q]* duh[q]; //squared norm of the gradient
323  * long double exponent=(p-2.0)/2*std::log(sqrGrad);
324  * computed_quantities[q](0) = std::sqrt(sqrGrad); // norm of the gradient
325  * computed_quantities[q](1)= std::exp(exponent); // multiplier
326  * }
327  * }
328  *
329  *
330  *
331  *
332  *
333  * template <int dim>
334  * std::vector<std::string>
335  * ComputeMultiplier<dim>::get_names() const
336  * {
337  * std::vector<std::string> solution_names;
338  * solution_names.push_back ("Gradient norm");
339  * solution_names.push_back ("Lagrange multiplier");
340  * return solution_names;
341  * }
342  *
343  *
344  * template <int dim>
345  * UpdateFlags
346  * ComputeMultiplier<dim>::get_needed_update_flags () const
347  * {
348  * return update_gradients;
349  * }
350  *
351  *
352  *
353  * template <int dim>
354  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
355  * ComputeMultiplier<dim>:: get_data_component_interpretation () const
356  * {
357  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
358  * interpretation;
359  * @endcode
360  *
361  * norm of the gradient
362  *
363  * @code
364  * interpretation.push_back (DataComponentInterpretation::component_is_scalar);
365  * @endcode
366  *
367  * Lagrange multiplier
368  *
369  * @code
370  * interpretation.push_back (DataComponentInterpretation::component_is_scalar);
371  * return interpretation;
372  * }
373  *
374  *
375  *
376  *
377  *
378  * @endcode
379  *
380  * ***************************************************************************************
381  *
382  * @code
383  * template <int dim>
384  * class ElastoplasticTorsion
385  * {
386  * public:
387  * ElastoplasticTorsion (ParameterHandler &);
388  * ~ElastoplasticTorsion ();
389  * void run ();
390  *
391  * private:
392  * void setup_system (const bool initial_step);
393  * void assemble_system ();
394  * bool solve (const int inner_it);
395  * void init_mesh ();
396  * void refine_mesh ();
397  * void set_boundary_values ();
398  * double phi (const double alpha) const;
399  * bool checkWolfe(double &alpha, double &phi_alpha) const;
400  * bool determine_step_length (const int inner_it);
401  * void print_it_message (const int counter, bool ks);
402  * void output_results (unsigned int refinement) const;
403  * void format_convergence_tables();
404  * void process_solution (const unsigned int cycle);
405  * void process_multiplier (const unsigned int cycle,const int iter,double time);
406  * double dual_error () const;
407  * double dual_infty_error () const;
408  * double W (double Du2) const;
409  * double Wp (double Du2) const;
410  * double G (double Du2) const;
411  *
412  *
413  *
414  * ParameterHandler &prm;
416  * DoFHandler<dim> dof_handler;
417  * AffineConstraints<double> hanging_node_constraints;
418  * SparsityPattern sparsity_pattern;
419  * SparseMatrix<double> system_matrix;
420  * ConvergenceTable convergence_table;
421  * ConvergenceTable dual_convergence_table;
422  * Vector<double> present_solution;
423  * Vector<double> newton_update;
424  * Vector<double> system_rhs;
425  * Vector<double> grad_norm;
427  *
428  *
429  * double step_length,phi_zero,phi_alpha,phip,phip_zero;
430  * double old_step,old_phi_zero,old_phip;
431  * double L2_error;
432  * double H1_error;
433  * double Linfty_error;
434  * double dual_L1_error;
435  * double dual_L_infty_error;
436  * FE_Q<dim> fe;
437  * double p;
438  * double line_search_tolerence; // c_1 in Nocedal & Wright
439  * unsigned int dir_id;
440  * std::string elements;
441  * std::string Method;
442  *
443  * };
444  *
445  * /*******************************************************************************************/
446  * @endcode
447  *
448  * Boundary condition
449  *
450 
451  *
452  *
453  * @code
454  * template <int dim>
455  * class BoundaryValues : public Function<dim>
456  * {
457  * public:
458  * BoundaryValues () : Function<dim>() {}
459  *
460  * virtual double value (const Point<dim> &p,
461  * const unsigned int component = 0) const;
462  * };
463  *
464  *
465  * template <int dim>
466  * double BoundaryValues<dim>::value (const Point<dim> &/*pto*/,
467  * const unsigned int /*component*/) const
468  * {
469  * @endcode
470  *
471  * could be anything else (theory works provided |Dg|_infty < 1/2)
472  *
473  * @code
474  * return 0.0;
475  *
476  * /* A challenging BC leading to overdetermined problems
477  * it is regulated by the parameter 0<eta<1.
478  * eta closer to 1 leads to more difficult problems.
479  *
480  * double pii=numbers::PI;
481  * double theta=std::atan2(p[1],p[0])+pii;
482  * double eta=0.9;
483  *
484  * if (theta <= 0.5)
485  * return eta*(theta*theta);
486  * else if ((theta >0.5) & (theta<= pii-0.5))
487  * return eta*(theta-0.25);
488  * else if ((theta>pii-0.5)&(theta<= pii+0.5))
489  * return eta*(pii-0.75-(theta-(pii-0.5))*(theta-(pii+0.5)));
490  * else if ((theta>pii+0.5)&(theta<= 2*pii-0.5))
491  * return eta*((2*pii-theta)-0.25);
492  * else
493  * return eta*((theta-2*pii)*(theta-2*pii) );*/
494  * }
495  *
496  *
497  *
498  * /******************************************************************************/
499  * @endcode
500  *
501  * Right-Hand Side
502  *
503  * @code
504  * template <int dim>
505  * class RightHandSide : public Function<dim>
506  * {
507  * public:
508  * RightHandSide () : Function<dim>() {}
509  * virtual double value (const Point<dim> &p,
510  * const unsigned int component = 0) const;
511  * };
512  *
513  * template <int dim>
514  * double RightHandSide<dim>::value (const Point<dim> &/*p*/,
515  * const unsigned int /*component*/) const
516  * {
517  * @endcode
518  *
519  * set to constant = 4, for which explicit solution to compare exists
520  * could be anything
521  *
522  * @code
523  * double return_value = 4.0;
524  * return return_value;
525  * }
526  *
527  *
528  *
529  * /*******************************************************************/
530  * @endcode
531  *
532  * The ElastoplasticTorsion class implementation
533  *
534 
535  *
536  * Constructor of the class
537  *
538  * @code
539  * template <int dim>
540  * ElastoplasticTorsion<dim>::ElastoplasticTorsion (ParameterHandler &param):
541  * prm(param),
542  * dof_handler (triangulation),
543  * L2_error(1.0),
544  * H1_error(1.0),
545  * Linfty_error(1.0),
546  * dual_L1_error(1.0),
547  * dual_L_infty_error(1.0),
548  * fe(2)
549  * {
550  * prm.enter_subsection ("Global Parameters");
551  * p=prm.get_double("p");
552  * prm.leave_subsection ();
553  * prm.enter_subsection ("Algorithm Parameters");
554  * line_search_tolerence=prm.get_double("line_search_tolerence");
555  * dir_id=prm.get_integer("Descent_direction");
556  * prm.leave_subsection ();
557  * if (fe.degree==1)
558  * elements="P1";
559  * else elements="P2";
560  *
561  * if (dir_id==0)
562  * Method="Precond";
563  * else
564  * Method="Newton";
565  * }
566  *
567  *
568  *
569  * template <int dim>
570  * ElastoplasticTorsion<dim>::~ElastoplasticTorsion ()
571  * {
572  * dof_handler.clear ();
573  * }
574  *
575  * /*****************************************************************************************/
576  * @endcode
577  *
578  * print iteration message
579  *
580 
581  *
582  *
583  * @code
584  * template <int dim>
585  * void ElastoplasticTorsion<dim>::print_it_message (const int counter, bool ks)
586  * {
587  * if (ks)
588  * {
589  * process_solution (counter);
590  * std::cout << "iteration="<< counter+1 << " J(u_h)= "<< phi_zero << ", H1 error: "
591  * << H1_error <<", W0-1,infty error: "<< Linfty_error<< " J'(u_h)(w)= "<< phip
592  * << ", |J'(u_h)|= "<< system_rhs.l2_norm()<<std::endl;
593  * }
594  * else
595  * {
596  * std::cout << "iteration= " << counter+1 << " J(u_h)= "
597  * << phi_alpha << " J'(u_h)= "<< phip<<std::endl;
598  * }
599  * }
600  *
601  *
602  * /*****************************************************************************************/
603  * @endcode
604  *
605  * Convergence Tables
606  *
607 
608  *
609  *
610 
611  *
612  *
613  * @code
614  * /*************************************************************/
615  * @endcode
616  *
617  * formating
618  *
619 
620  *
621  *
622  * @code
623  * template <int dim>
624  * void ElastoplasticTorsion<dim>::format_convergence_tables()
625  * {
626  * convergence_table.set_precision("L2", 3);
627  * convergence_table.set_precision("H1", 3);
628  * convergence_table.set_precision("Linfty", 3);
629  * convergence_table.set_precision("function value", 3);
630  * convergence_table.set_precision("derivative", 3);
631  * dual_convergence_table.set_precision("dual_L1", 3);
632  * dual_convergence_table.set_precision("dual_Linfty", 3);
633  * dual_convergence_table.set_precision("L2", 3);
634  * dual_convergence_table.set_precision("H1", 3);
635  * dual_convergence_table.set_precision("Linfty", 3);
636  * convergence_table.set_scientific("L2", true);
637  * convergence_table.set_scientific("H1", true);
638  * convergence_table.set_scientific("Linfty", true);
639  * convergence_table.set_scientific("function value", true);
640  * convergence_table.set_scientific("derivative", true);
641  * dual_convergence_table.set_scientific("dual_L1", true);
642  * dual_convergence_table.set_scientific("dual_Linfty", true);
643  * dual_convergence_table.set_scientific("L2", true);
644  * dual_convergence_table.set_scientific("H1", true);
645  * dual_convergence_table.set_scientific("Linfty", true);
646  *
647  * }
648  *
649  * /****************************************/
650  * @endcode
651  *
652  * fill-in entry for the solution
653  *
654  * @code
655  * template <int dim>
656  * void ElastoplasticTorsion<dim>::process_solution (const unsigned int it)
657  * {
658  * Vector<float> difference_per_cell (triangulation.n_active_cells());
659  *
660  * @endcode
661  *
662  * compute L2 error (save to difference_per_cell)
663  *
664  * @code
665  * VectorTools::integrate_difference (dof_handler,present_solution,
666  * Solution<dim>(),difference_per_cell,QGauss<dim>(3),VectorTools::L2_norm);
667  * L2_error = difference_per_cell.l2_norm();
668  *
669  * @endcode
670  *
671  * compute H1 error (save to difference_per_cell)
672  *
673  * @code
674  * VectorTools::integrate_difference (dof_handler,present_solution,Solution<dim>(),
675  * difference_per_cell,QGauss<dim>(3),VectorTools::H1_seminorm);
676  * H1_error = difference_per_cell.l2_norm();
677  *
678  * @endcode
679  *
680  * compute W1infty error (save to difference_per_cell)
681  *
682  * @code
683  * const QTrapez<1> q_trapez;
684  * const QIterated<dim> q_iterated (q_trapez, 5);
685  * VectorTools::integrate_difference (dof_handler,present_solution,Solution<dim>(),
686  * difference_per_cell,q_iterated,VectorTools::W1infty_seminorm);
687  * Linfty_error = difference_per_cell.linfty_norm();
688  *
689  *
690  * convergence_table.add_value("cycle", it);
691  * convergence_table.add_value("p", p);
692  * convergence_table.add_value("L2", L2_error);
693  * convergence_table.add_value("H1", H1_error);
694  * convergence_table.add_value("Linfty", Linfty_error);
695  * convergence_table.add_value("function value", phi_alpha);
696  * convergence_table.add_value("derivative", phip);
697  * }
698  *
699  *
700  * /***************************************/
701  * @endcode
702  *
703  * fill-in entry for the multiplier
704  *
705  * @code
706  * template <int dim>
707  * void ElastoplasticTorsion<dim>::process_multiplier (const unsigned int cycle, const int iter,double time)
708  * {
709  * const unsigned int n_active_cells=triangulation.n_active_cells();
710  * const unsigned int n_dofs=dof_handler.n_dofs();
711  * dual_L1_error=dual_error();
712  * dual_L_infty_error=dual_infty_error();
713  *
714  *
715  * dual_convergence_table.add_value("cycle", cycle);
716  * dual_convergence_table.add_value("p", p);
717  * dual_convergence_table.add_value("iteration_number", iter);
718  * dual_convergence_table.add_value("cpu_time", time);
719  * dual_convergence_table.add_value("cells", n_active_cells);
720  * dual_convergence_table.add_value("dofs", n_dofs);
721  * dual_convergence_table.add_value("L2", L2_error);
722  * dual_convergence_table.add_value("H1", H1_error);
723  * dual_convergence_table.add_value("Linfty", Linfty_error);
724  * dual_convergence_table.add_value("dual_L1", dual_L1_error);
725  * dual_convergence_table.add_value("dual_Linfty", dual_L_infty_error);
726  *
727  * }
728  *
729  *
730  *
731  *
732  * /****************************************************************************************/
733  * @endcode
734  *
735  * ElastoplasticTorsion::setup_system
736  * unchanged from @ref step_15 "step-15"
737  *
738 
739  *
740  *
741  * @code
742  * template <int dim>
743  * void ElastoplasticTorsion<dim>::setup_system (const bool initial_step)
744  * {
745  * if (initial_step)
746  * {
747  * dof_handler.distribute_dofs (fe);
748  * present_solution.reinit (dof_handler.n_dofs());
749  * grad_norm.reinit (dof_handler.n_dofs());
750  * lambda.reinit (dof_handler.n_dofs());
751  *
752  * hanging_node_constraints.clear ();
754  * hanging_node_constraints);
755  * hanging_node_constraints.close ();
756  * }
757  *
758  *
759  * @endcode
760  *
761  * The remaining parts of the function
762  *
763 
764  *
765  *
766  * @code
767  * newton_update.reinit (dof_handler.n_dofs());
768  * system_rhs.reinit (dof_handler.n_dofs());
769  * DynamicSparsityPattern c_sparsity(dof_handler.n_dofs());
770  * DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
771  * hanging_node_constraints.condense (c_sparsity);
772  * sparsity_pattern.copy_from(c_sparsity);
773  * system_matrix.reinit (sparsity_pattern);
774  * }
775  *
776  * /***************************************************************************************/
777  * /* the coeffcients W, W' and G defining the problem.
778  *
779  * Min_u \int W(|Du|^2) dx
780  *
781  * They must be consistent as G(s)=W'(s)+2s W''(s) for any s>0.
782  * recall that they receive the SQUARED gradient. */
783  *
784  * template <int dim>
785  * double ElastoplasticTorsion<dim>::W (double Du2) const
786  * {
787  * return Du2;
788  * }
789  *
790  * template <int dim>
791  * double ElastoplasticTorsion<dim>::Wp (double /*Du2*/) const
792  * {
793  * return 1.0;
794  * }
795  *
796  * template <int dim>
797  * double ElastoplasticTorsion<dim>::G (double /*Du2*/) const
798  * {
799  * return 1.0;
800  * }
801  * /***************************************************************************************/
802  *
803  * template <int dim>
804  * void ElastoplasticTorsion<dim>::assemble_system ()
805  * {
806  *
807  *
808  * const QGauss<dim> quadrature_formula(3);
809  * const RightHandSide<dim> right_hand_side;
810  * system_matrix = 0;
811  * system_rhs = 0;
812  *
813  * FEValues<dim> fe_values (fe, quadrature_formula,
814  * update_gradients |
815  * update_values |
818  *
819  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
820  * const unsigned int n_q_points = quadrature_formula.size();
821  *
822  * FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
823  * Vector<double> cell_rhs (dofs_per_cell);
824  *
825  * std::vector<Tensor<1, dim> > old_solution_gradients(n_q_points);
826  * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
827  *
828  *
830  * cell = dof_handler.begin_active(),
831  * endc = dof_handler.end();
832  * for (; cell!=endc; ++cell)
833  * {
834  * cell_matrix = 0;
835  * cell_rhs = 0;
836  *
837  * fe_values.reinit (cell);
838  * fe_values.get_function_gradients(present_solution,
839  * old_solution_gradients);
840  *
841  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
842  * {
843  * long double coeff=0.0;
844  * long double a=old_solution_gradients[q_point] * old_solution_gradients[q_point];
845  * long double exponent=(p-2.0)/2*std::log(a);
846  * coeff= std::exp( exponent);
847  * for (unsigned int i=0; i<dofs_per_cell; ++i)
848  * {
849  * for (unsigned int j=0; j<dofs_per_cell; ++j)
850  * {
851  * if (dir_id==1)
852  * {
853  * cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * fe_values.shape_grad(j, q_point)
854  * * (G(a)+(p-1.0)*coeff) * fe_values.JxW(q_point);
855  * }
856  * else
857  * {
858  * cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * fe_values.shape_grad(j, q_point)
859  * * (Wp(a)+coeff)
860  * * fe_values.JxW(q_point);
861  * }
862  * }
863  *
864  * cell_rhs(i) -= ( fe_values.shape_grad(i, q_point)
865  * * old_solution_gradients[q_point]
866  * * (Wp(a)+coeff)
867  * -right_hand_side.value(fe_values.quadrature_point(q_point))
868  * *fe_values.shape_value(i, q_point)
869  * )
870  * * fe_values.JxW(q_point);
871  * }
872  * }
873  *
874  * cell->get_dof_indices (local_dof_indices);
875  * for (unsigned int i=0; i<dofs_per_cell; ++i)
876  * {
877  * for (unsigned int j=0; j<dofs_per_cell; ++j)
878  * system_matrix.add (local_dof_indices[i],
879  * local_dof_indices[j],
880  * cell_matrix(i,j));
881  *
882  * system_rhs(local_dof_indices[i]) += cell_rhs(i);
883  * }
884  * }
885  *
886  * hanging_node_constraints.condense (system_matrix);
887  * hanging_node_constraints.condense (system_rhs);
888  *
889  * std::map<types::global_dof_index,double> boundary_values;
891  * 0,
892  * ZeroFunction<dim>(),
893  * boundary_values);
894  * MatrixTools::apply_boundary_values (boundary_values,
895  * system_matrix,
896  * newton_update,
897  * system_rhs);
898  * }
899  *
900  *
901  *
902  *
903  * /********************************** Refine Mesh ****************************************/
904  * @endcode
905  *
906  * unchanged from @ref step_15 "step-15"
907  *
908 
909  *
910  *
911  * @code
912  * template <int dim>
913  * void ElastoplasticTorsion<dim>::refine_mesh ()
914  * {
915  * using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
916  *
917  * Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
918  * KellyErrorEstimator<dim>::estimate (dof_handler,
919  * QGauss<dim-1>(3),
920  * FunctionMap(),
921  * present_solution,
922  * estimated_error_per_cell);
923  *
924  * prm.enter_subsection ("Mesh & Refinement Parameters");
925  * const double top_fraction=prm.get_double("top_fraction_of_cells");
926  * const double bottom_fraction=prm.get_double("bottom_fraction_of_cells");
927  * prm.leave_subsection ();
929  * estimated_error_per_cell,
930  * top_fraction, bottom_fraction);
931  *
932  * triangulation.prepare_coarsening_and_refinement ();
933  * SolutionTransfer<dim> solution_transfer(dof_handler);
934  * solution_transfer.prepare_for_coarsening_and_refinement(present_solution);
935  * triangulation.execute_coarsening_and_refinement();
936  * dof_handler.distribute_dofs(fe);
937  * Vector<double> tmp(dof_handler.n_dofs());
938  * solution_transfer.interpolate(present_solution, tmp);
939  * present_solution = tmp;
940  * set_boundary_values ();
941  * hanging_node_constraints.clear();
942  *
944  * hanging_node_constraints);
945  * hanging_node_constraints.close();
946  * hanging_node_constraints.distribute (present_solution);
947  * setup_system (false);
948  * }
949  *
950  *
951  * /*******************************************************************************************/
952  * @endcode
953  *
954  * Dump the norm of the gradient and the lagrange multiplier in vtu format for visualization
955  *
956  * @code
957  * template <int dim>
958  * void ElastoplasticTorsion<dim>::output_results (unsigned int counter) const
959  * {
960  * @endcode
961  *
962  * multiplier object contains both |Du| and lambda.
963  *
964  * @code
965  * ComputeMultiplier<dim> multiplier(p);
966  * DataOut<dim> data_out;
967  *
968  * data_out.attach_dof_handler (dof_handler);
969  * data_out.add_data_vector (present_solution, "solution");
970  * data_out.add_data_vector (present_solution, multiplier);
971  * data_out.build_patches ();
972  * std::ostringstream p_str;
973  * p_str << p<<"-cycle-"<<counter;
974  * std::string str = p_str.str();
975  * const std::string filename = "solution-" + str+".vtu";
976  * std::ofstream output (filename.c_str());
977  * data_out.write_vtu (output);
978  * }
979  *
980  * /********************************************************************************************/
981  * @endcode
982  *
983  * unchanged from @ref step_15 "step-15"
984  *
985  * @code
986  * template <int dim>
987  * void ElastoplasticTorsion<dim>::set_boundary_values ()
988  * {
989  * std::map<types::global_dof_index, double> boundary_values;
991  * 0,
992  * BoundaryValues<dim>(),
993  * boundary_values);
994  * for (std::map<types::global_dof_index, double>::const_iterator
995  * bp = boundary_values.begin();
996  * bp != boundary_values.end(); ++bp)
997  * present_solution(bp->first) = bp->second;
998  * }
999  *
1000  *
1001  * /****************************************************************************************/
1002  * @endcode
1003  *
1004  * COMPUTE @f$\phi(\alpha)=J_p(u_h+\alpha w)@f$
1005  *
1006  * @code
1007  * template <int dim>
1008  * double ElastoplasticTorsion<dim>::phi (const double alpha) const
1009  * {
1010  * double obj = 0.0;
1011  * const RightHandSide<dim> right_hand_side;
1012  * Vector<double> evaluation_point (dof_handler.n_dofs());
1013  * evaluation_point = present_solution; // copy of u_h
1014  * evaluation_point.add (alpha, newton_update); // u_{n+1}=u_n+alpha w_n
1015  *
1016  * const QGauss<dim> quadrature_formula(3);
1017  * FEValues<dim> fe_values (fe, quadrature_formula,
1018  * update_gradients |
1019  * update_values |
1021  * update_JxW_values);
1022  *
1023  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
1024  * const unsigned int n_q_points = quadrature_formula.size();
1025  *
1026  * Vector<double> cell_residual (dofs_per_cell);
1027  * std::vector<Tensor<1, dim> > gradients(n_q_points);
1028  * std::vector<double> values(n_q_points);
1029  *
1030  *
1031  * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1032  *
1034  * cell = dof_handler.begin_active(),
1035  * endc = dof_handler.end();
1036  * for (; cell!=endc; ++cell)
1037  * {
1038  * cell_residual = 0;
1039  * fe_values.reinit (cell);
1040  * fe_values.get_function_gradients (evaluation_point, gradients);
1041  * fe_values.get_function_values (evaluation_point, values);
1042  *
1043  *
1044  * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
1045  * {
1046  * double Du2=gradients[q_point] * gradients[q_point]; // Du2=|Du|^2
1047  * double penalty;
1048  * if (Du2<1.0e-10)
1049  * penalty=0.0;
1050  * else
1051  * penalty=std::pow(Du2,p/2.0); // penalty=|Du|^p
1052  *
1053  * @endcode
1054  *
1055  * obj+= 1/2 W(|Du|^2)+1/p |Du|^p -fu (see (1))
1056  *
1057  * @code
1058  * obj+=(
1059  * (0.5*W(Du2)+penalty/p)- right_hand_side.value(fe_values.quadrature_point(q_point))*values[q_point]
1060  * ) * fe_values.JxW(q_point);
1061  * }
1062  *
1063  * }
1064  *
1065  * return obj;
1066  * }
1067  *
1068  *
1069  * /***************************************************************************************************/
1070  * @endcode
1071  *
1072  * Compute L^1 error norm of Lagrange Multiplier
1073  * with respect to exact solution (cf. Alvarez & Flores, 2015)
1074  *
1075 
1076  *
1077  *
1078  * @code
1079  * template <int dim>
1080  * double ElastoplasticTorsion<dim>::dual_error () const
1081  * {
1082  * double obj = 0.0;
1083  *
1084  * const QGauss<dim> quadrature_formula(3);
1085  * FEValues<dim> fe_values (fe, quadrature_formula,
1086  * update_gradients |
1088  * update_JxW_values);
1089  *
1090  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
1091  * const unsigned int n_q_points = quadrature_formula.size();
1092  *
1093  * Vector<double> cell_residual (dofs_per_cell);
1094  * std::vector<Tensor<1, dim> > gradients(n_q_points);
1095  *
1096  * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1097  *
1099  * cell = dof_handler.begin_active(),
1100  * endc = dof_handler.end();
1101  * for (; cell!=endc; ++cell)
1102  * {
1103  * cell_residual = 0;
1104  * fe_values.reinit (cell);
1105  * fe_values.get_function_gradients (present_solution, gradients);
1106  *
1107  * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
1108  * {
1109  * double coeff=gradients[q_point] * gradients[q_point] ;
1110  * if (coeff<1.0e-15)
1111  * coeff=0.0;
1112  * else
1113  * coeff=std::pow(coeff,(p-2.0)/2.0); // |Du_p|^(p-2)
1114  *
1115  * double r=std::sqrt(fe_values.quadrature_point(q_point).square());
1116  * double exact=0;
1117  * if (r>0.5)
1118  * exact= 2*r-1;
1119  *
1120  * obj+=( std::abs(coeff-exact) ) * fe_values.JxW(q_point);
1121  * }
1122  *
1123  * }
1124  *
1125  * return obj;
1126  * }
1127  *
1128  * /*******************************************************************************************/
1129  * @endcode
1130  *
1131  * Compute L^infinity error norm of Lagrange Multiplier
1132  * with respect to exact solution (cf. Alvarez & Flores, 2015)
1133  *
1134 
1135  *
1136  *
1137  * @code
1138  * template <int dim>
1139  * double ElastoplasticTorsion<dim>::dual_infty_error () const
1140  * {
1141  * double obj = 0.0;
1142  * const QTrapez<1> q_trapez;
1143  * const QIterated<dim> quadrature_formula (q_trapez, 10);
1144  *
1145  * FEValues<dim> fe_values (fe, quadrature_formula,
1146  * update_gradients |
1148  *
1149  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
1150  * const unsigned int n_q_points = quadrature_formula.size();
1151  *
1152  * Vector<double> cell_residual (dofs_per_cell);
1153  * std::vector<Tensor<1, dim> > gradients(n_q_points);
1154  *
1155  * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1156  *
1158  * cell = dof_handler.begin_active(),
1159  * endc = dof_handler.end();
1160  * for (; cell!=endc; ++cell)
1161  * {
1162  * cell_residual = 0;
1163  * fe_values.reinit (cell);
1164  * fe_values.get_function_gradients (present_solution, gradients);
1165  *
1166  * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
1167  * {
1168  * long double sqdGrad=gradients[q_point] * gradients[q_point] ;
1169  * double r=std::sqrt(fe_values.quadrature_point(q_point).square());
1170  * double exact=0;
1171  * if (r>0.5)
1172  * exact= 2*r-1.0;
1173  * @endcode
1174  *
1175  * compute |Du|^(p-2) as exp(p-2/2*log(Du^2))
1176  *
1177  * @code
1178  * long double exponent=(p-2.0)/2*std::log(sqdGrad);
1179  * long double coeff=std::exp(exponent);
1180  *
1181  * if (std::abs(coeff-exact)>obj )
1182  * obj=std::abs(coeff-exact);
1183  * }
1184  *
1185  * }
1186  *
1187  * return obj;
1188  * }
1189  *
1190  * /*****************************************************************************************/
1191  * @endcode
1192  *
1193  * check whether putative @ref step_length "step-length" satisfies sufficient decrease conditions
1194  *
1195  * @code
1196  * template <int dim>
1197  * bool ElastoplasticTorsion<dim>::checkWolfe(double &alpha, double &phi_alpha) const
1198  * {
1199  * if (phi_alpha< phi_zero+line_search_tolerence*phip*alpha )
1200  * return true;
1201  * else
1202  * return false;
1203  * }
1204  *
1205  *
1206  * /*****************************************************************************************/
1207  * @endcode
1208  *
1209  * Find a @ref step_length "step-length" satisfying sufficient decrease condition by line-search
1210  * uses quadratic interpolation
1211  *
1212 
1213  *
1214  *
1215  * @code
1216  * template <int dim>
1217  * bool ElastoplasticTorsion<dim>::determine_step_length(const int inner_it)
1218  * {
1219  * unsigned int it=0;
1220  * bool done;
1221  * double alpha,nalpha;
1222  * prm.enter_subsection ("Algorithm Parameters");
1223  * const unsigned int max_LS_it=prm.get_integer("max_LS_it");
1224  * double init_SL=prm.get_double("init_step_length");
1225  * prm.leave_subsection ();
1226  * if (inner_it==0)
1227  * alpha=init_SL;
1228  * else
1229  * {
1230  * alpha=std::min(1.45*old_step*old_phip/phip,1.0);
1231  * }
1232  * phi_alpha=phi(alpha);
1233  * std::cerr << "Step length=" << alpha << ", Value= " << phi_alpha;
1234  * @endcode
1235  *
1236  * check if @ref step_size "step-size" satisfies sufficient decrease condition
1237  *
1238  * @code
1239  * done=checkWolfe(alpha,phi_alpha);
1240  * if (done)
1241  * std::cerr << " accepted" << std::endl;
1242  * else
1243  * std::cerr << " rejected" ;
1244  *
1245  * while ((!done) & (it<max_LS_it))
1246  * {
1247  * @endcode
1248  *
1249  * new try obtained by quadratic interpolation
1250  *
1251  * @code
1252  * nalpha=-(phip*alpha*alpha)/(2*(phi_alpha-phi_zero-phip*alpha));
1253  *
1254  * if (nalpha<1e-3*alpha || std::abs(nalpha-alpha)/alpha<1e-8)
1255  * nalpha=alpha/2;
1256  * else if ( phi_alpha-phi_zero>1e3*std::abs(phi_zero) )
1257  * nalpha=alpha/10;
1258  * alpha=nalpha;
1259  * phi_alpha=phi(alpha);
1260  * done=checkWolfe(alpha,phi_alpha);
1261  * if (done)
1262  * std::cerr << ", finished with steplength= "<< alpha<< ", fcn value= "<< phi_alpha<<std::endl;
1263  * it=it+1;
1264  * }
1265  * if (!done)
1266  * {
1267  * std::cerr << ", max. no. of iterations reached wiht steplength= "<< alpha
1268  * << ", fcn value= "<< phi_alpha<<std::endl;
1269  * return false;
1270  * }
1271  * else
1272  * {
1273  * step_length=alpha;
1274  * return true;
1275  * }
1276  *
1277  * }
1278  *
1279  * /**************************************************************************************************/
1280  * @endcode
1281  *
1282  * ElastoplasticTorsion::init_mesh()
1283  *
1284 
1285  *
1286  *
1287  * @code
1288  * template <int dim>
1289  * void ElastoplasticTorsion<dim>::init_mesh ()
1290  * {
1291  * @endcode
1292  *
1293  * get parameters
1294  *
1295  * @code
1296  * prm.enter_subsection ("Mesh & Refinement Parameters");
1297  * const int domain_id=prm.get_integer("Code for the domain");
1298  * const int init_ref=prm.get_integer("No of initial refinements");
1299  * prm.leave_subsection ();
1300  *
1301  *
1302  * if (domain_id==0)
1303  * {
1304  * @endcode
1305  *
1306  * For the unit disk around the origin
1307  *
1308  * @code
1310  * static const SphericalManifold<dim> boundary;
1311  * triangulation.set_manifold (0, boundary);
1312  * }
1313  * else if (domain_id==1)
1314  * {
1315  * @endcode
1316  *
1317  * For the unit square
1318  *
1319  * @code
1321  * }
1322  * else if (domain_id==2)
1323  * {
1324  * /* For Glowinski's domain
1325  * ___ ___ __ 1
1326  * | |__| | __ .8
1327  * | |
1328  * | |
1329  * |__________| __ 0
1330  *
1331  * | | | |
1332  * 0 .4 .6 1
1333  *
1334  * */
1335  * Triangulation<dim> tria1;
1336  * Triangulation<dim> tria2;
1337  * Triangulation<dim> tria3;
1338  * Triangulation<dim> tria4;
1339  * Triangulation<dim> tria5;
1340  * Triangulation<dim> tria6;
1341  * GridGenerator::hyper_rectangle(tria1, Point<2>(0.0,0.0), Point<2>(0.4,0.8));
1342  * GridGenerator::hyper_rectangle(tria2, Point<2>(0.0,0.8), Point<2>(0.4,1.0));
1343  * GridGenerator::hyper_rectangle(tria3, Point<2>(0.4,0.0), Point<2>(0.6,0.8));
1344  * GridGenerator::hyper_rectangle(tria4, Point<2>(0.6,0.0), Point<2>(1.0,0.8));
1345  * GridGenerator::hyper_rectangle(tria5, Point<2>(0.6,0.8), Point<2>(1.0,1.0));
1346  * GridGenerator::merge_triangulations (tria1, tria2, tria6);
1347  * GridGenerator::merge_triangulations (tria6, tria3, tria6);
1348  * GridGenerator::merge_triangulations (tria6, tria4, tria6);
1350  * }
1351  * @endcode
1352  *
1353  * perform initial refinements
1354  *
1355  * @code
1356  * triangulation.refine_global(init_ref);
1357  * }
1358  *
1359  * /**************************************************************************************************/
1360  * @endcode
1361  *
1362  * ElastoplasticTorsion::solve(inner_it)
1363  * Performs one inner iteration
1364  *
1365 
1366  *
1367  *
1368  * @code
1369  * template <int dim>
1370  * bool ElastoplasticTorsion<dim>::solve (const int inner_it)
1371  * {
1372  * prm.enter_subsection ("Algorithm Parameters");
1373  * const unsigned int max_CG_it=prm.get_integer("Max_CG_it");
1374  * const double CG_tol=prm.get_double("CG_tol");
1375  * prm.leave_subsection ();
1376  *
1377  * SolverControl solver_control (max_CG_it,CG_tol);
1378  * SolverCG<> solver (solver_control);
1379  *
1380  * PreconditionSSOR<> preconditioner;
1381  * preconditioner.initialize(system_matrix,0.25);
1382  *
1383  * solver.solve (system_matrix, newton_update, system_rhs,
1384  * preconditioner);
1385  * hanging_node_constraints.distribute (newton_update);
1386  * /****** save current quantities for line-search **** */
1387  * @endcode
1388  *
1389  * Recall that phi(alpha)=J(u+alpha w)
1390  *
1391  * @code
1392  * old_step=step_length;
1393  * old_phi_zero=phi_zero;
1394  * phi_zero=phi(0); // phi(0)=J(u)
1395  * old_phip=phip;
1396  * phip=-1.0*(newton_update*system_rhs); //phi'(0)=J'(u) *w, rhs=-J'(u).
1397  * if (inner_it==0)
1398  * phip_zero=phip;
1399  *
1400  * if (phip>0) // this should not happen, step back
1401  * {
1402  * std::cout << "Not a descent direction!" <<std::endl;
1403  * present_solution.add (-1.0*step_length, newton_update);
1404  * step_length=step_length/2;
1405  * phip=old_phip;
1406  * return false;
1407  * }
1408  * else
1409  * {
1410  * if (determine_step_length(inner_it))
1411  * {
1412  * @endcode
1413  *
1414  * update u_{n+1}=u_n+alpha w_n
1415  *
1416  * @code
1417  * present_solution.add (step_length, newton_update);
1418  * return true;
1419  * }
1420  * else return false;
1421  * }
1422  * }
1423  *
1424  *
1425  *
1426  * /*************************************************************************************************************/
1427  * @endcode
1428  *
1430  *
1431  * @code
1432  * template <int dim>
1434  * {
1435  *
1436  * @endcode
1437  *
1438  * get parameters
1439  *
1440  * @code
1441  * prm.enter_subsection ("Mesh & Refinement Parameters");
1442  * const int adapt_ref=prm.get_integer("No of adaptive refinements");
1443  * prm.leave_subsection ();
1444  * prm.enter_subsection ("Algorithm Parameters");
1445  * const int max_inner=prm.get_integer("Max_inner");
1446  * const double eps=prm.get_double("eps");
1447  * const double hi_eps=prm.get_double("hi_eps");
1448  * const int hi_th=prm.get_integer("hi_th");
1449  * const double init_p=prm.get_double("init_p");
1450  * const double delta_p=prm.get_double("delta_p");
1451  * prm.leave_subsection ();
1452  * prm.enter_subsection ("Global Parameters");
1453  * bool known_solution=prm.get_bool("known_solution");
1454  * double actual_p=prm.get_double("p");
1455  * prm.leave_subsection ();
1456  * /************************/
1457  *
1458  * @endcode
1459  *
1460  * init Timer
1461  *
1462  * @code
1463  * Timer timer;
1464  * double ptime=0.0;
1465  * timer.start ();
1466  *
1467  * @endcode
1468  *
1469  * initalize mesh for the selected domain
1470  *
1471  * @code
1472  * init_mesh();
1473  *
1474  * @endcode
1475  *
1476  * setup FE space
1477  *
1478  * @code
1479  * setup_system (true);
1480  * set_boundary_values ();
1481  *
1482  * @endcode
1483  *
1484  * init counters
1485  *
1486  * @code
1487  * int global_it=0; // Total inner iterations (counting both loops)
1488  * int cycle=0; // Total outer iterations (counting both loops)
1489  * int refinement = 0; // Refinements performed (adaptive) = outer iterations 2nd loop
1490  *
1491  *
1492  * @endcode
1493  *
1494  * prepare to start first loop
1495  *
1496  * @code
1497  * p=init_p;
1498  * bool well_solved=true;
1499  *
1500  * /***************************** First loop ***********************************/
1501  * /****************** Prepare initial condition using increasing p *************************/
1502  * while (p<actual_p) // outer iteration, increasing p.
1503  * {
1504  * std::cout <<"--Preparing initial condition with p="<<p<<" iter.= " << global_it<< " .-- "<< std::endl;
1505  * timer.restart();
1506  * for (int inner_iteration=0; inner_iteration<max_inner; ++inner_iteration,++global_it)
1507  * {
1508  * assemble_system ();
1509  * well_solved=solve (inner_iteration);
1510  * print_it_message (global_it, known_solution);
1511  * if (
1512  * ((system_rhs.l2_norm()/std::sqrt(system_rhs.size()) <1e-4) & (cycle<1)) |
1513  * ((system_rhs.l2_norm()/std::sqrt(system_rhs.size()) <1e-5) & (cycle>=1)) |
1514  * !well_solved
1515  * )
1516  * break;
1517  * }
1518  * ptime=timer.cpu_time();
1519  * if (well_solved)
1520  * output_results (cycle);
1521  *
1522  * if (known_solution)
1523  * {
1524  * process_multiplier(cycle,global_it,ptime);
1525  * @endcode
1526  *
1527  * dual_convergence_table.write_tex(dual_error_table_file);
1528  *
1529  * @code
1530  * }
1531  * refine_mesh();
1532  * cycle++;
1533  * p+=delta_p;
1534  * }
1535  * /*************************** first loop finished ********************/
1536  *
1537  *
1538  * @endcode
1539  *
1540  * prepare for second loop
1541  *
1542  * @code
1543  * p=actual_p;
1544  * well_solved=true;
1545  *
1546  *
1547  * /***************************** Second loop *********************************/
1548  * /**************************** Solve problem for target p *********************************/
1549  *
1550  * std::cout << "============ Solving problem with p=" <<p << " ==================" << std::endl;
1551  * /***** Outer iteration - refining mesh ******************/
1552  * while ((cycle<adapt_ref) & well_solved)
1553  * {
1554  * timer.restart();
1555  * @endcode
1556  *
1557  * inner iteration
1558  *
1559  * @code
1560  * for (int inner_iteration=0; inner_iteration<max_inner; ++inner_iteration,++global_it)
1561  * {
1562  * assemble_system ();
1563  * well_solved=solve (inner_iteration);
1564  * print_it_message (global_it, known_solution);
1565  *
1566  * if (
1567  * ((system_rhs.l2_norm()/std::sqrt(system_rhs.size()) < eps) & (refinement<hi_th)) |
1568  * (( system_rhs.l2_norm()/ std::sqrt (system_rhs.size()) <hi_eps) | (!well_solved))
1569  * )
1570  * break;
1571  * }
1572  * @endcode
1573  *
1574  * inner iterations finished
1575  *
1576  * @code
1577  * ptime=timer.cpu_time();
1578  * if (well_solved)
1579  * output_results (cycle);
1580  *
1581  * @endcode
1582  *
1583  * compute and display error, if the explicit solution is known
1584  *
1585  * @code
1586  * if (known_solution)
1587  * {
1588  * process_multiplier(cycle,global_it,ptime);
1589  * std::cout << "finished with H1 error: " << H1_error << ", dual error (L1): "
1590  * << dual_L1_error << "dual error (L infty): "<<dual_L_infty_error <<std::endl;
1591  * }
1592  *
1593  * @endcode
1594  *
1595  * update counters
1596  *
1597  * @code
1598  * ++refinement;
1599  * ++cycle;
1600  * @endcode
1601  *
1602  * refine mesh
1603  *
1604  * @code
1605  * std::cout << "******** Refined mesh " << cycle << " ********" << std::endl;
1606  * refine_mesh();
1607  * }// second loop
1608  *
1609  * @endcode
1610  *
1611  * write convergence tables to file
1612  *
1613  * @code
1614  * if (known_solution)
1615  * {
1616  * format_convergence_tables();
1617  * std::string error_filename = "error"+Method+elements+".tex";
1618  * std::ofstream error_table_file(error_filename.c_str());
1619  * std::string dual_error_filename = "dual_error"+Method+elements+".tex";
1620  * std::ofstream dual_error_table_file(dual_error_filename.c_str());
1621  * convergence_table.write_tex(error_table_file);
1622  * dual_convergence_table.write_tex(dual_error_table_file);
1623  * }
1624  * }//run()
1625  *
1626  * }//namespace
1627  *
1628  * /**********************************************************************************************/
1629  * @endcode
1630  *
1631  * The main function
1632  *
1633  * @code
1634  * int main ()
1635  * {
1636  * try
1637  * {
1638  * using namespace dealii;
1639  * using namespace nsp;
1640  * deallog.depth_console (0);
1641  *
1642  * ParameterHandler prm;
1643  * ParameterReader param(prm);
1644  * param.read_parameters("EPT.prm");
1645  * ElastoplasticTorsion<2> ElastoplasticTorsionProblem(prm);
1646  * ElastoplasticTorsionProblem .run ();
1647  * }
1648  * catch (std::exception &exc)
1649  * {
1650  * std::cerr << std::endl << std::endl
1651  * << "----------------------------------------------------" << std::endl;
1652  * std::cerr << "Exception on processing: " << std::endl
1653  * << exc.what() << std::endl
1654  * << "Aborting!" << std::endl
1655  * << "----------------------------------------------------"
1656  * << std::endl;
1657  *
1658  * return 1;
1659  * }
1660  * catch (...)
1661  * {
1662  * std::cerr << std::endl << std::endl
1663  * << "----------------------------------------------------"
1664  * << std::endl;
1665  * std::cerr << "Unknown exception!" << std::endl
1666  * << "Aborting!" << std::endl
1667  * << "----------------------------------------------------"
1668  * << std::endl;
1669  * return 1;
1670  * }
1671  * return 0;
1672  * }
1673  * @endcode
1674 
1675 
1676 */
std::exp
inline ::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5368
DataPostprocessor
Definition: data_postprocessor.h:502
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
DataOutInterface::write_vtu
void write_vtu(std::ostream &out) const
Definition: data_out_base.cc:6864
VectorTools::L2_norm
@ L2_norm
Definition: vector_tools_common.h:113
SolverCG
Definition: solver_cg.h:98
FE_Q
Definition: fe_q.h:554
dealii
Definition: namespace_dealii.h:25
DataComponentInterpretation::component_is_scalar
@ component_is_scalar
Definition: data_component_interpretation.h:55
Differentiation::SD::OptimizerType::lambda
@ lambda
Triangulation< dim >
VectorTools::W1infty_seminorm
@ W1infty_seminorm
Definition: vector_tools_common.h:258
GridRefinement::refine_and_coarsen_fixed_number
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())
Definition: grid_refinement.cc:189
SparseMatrix< double >
DataOutBase::vtu
@ vtu
Definition: data_out_base.h:1610
Patterns::Bool
Definition: patterns.h:984
LAPACKSupport::L
static const char L
Definition: lapack_support.h:171
MatrixTools::apply_boundary_values
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
Definition: matrix_tools.cc:81
SphericalManifold
Definition: manifold_lib.h:231
VectorTools::integrate_difference
void integrate_difference(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const InVector &fe_function, const Function< spacedim, typename InVector::value_type > &exact_solution, OutVector &difference, const Quadrature< dim > &q, const NormType &norm, const Function< spacedim, double > *weight=nullptr, const double exponent=2.)
GridGenerator::hyper_rectangle
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Function::gradient
virtual Tensor< 1, dim, RangeNumberType > gradient(const Point< dim > &p, const unsigned int component=0) const
VectorizedArray::log
VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5390
LocalIntegrators::Advection::cell_matrix
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.)
Definition: advection.h:80
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
second
Point< 2 > second
Definition: grid_out.cc:4353
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1071
DoFHandler< dim >
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
deallog
LogStream deallog
Definition: logstream.cc:37
QIterated
Definition: quadrature.h:369
LocalIntegrators::Advection::cell_residual
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim >> &input, const ArrayView< const std::vector< double >> &velocity, double factor=1.)
Definition: advection.h:136
FunctionMap
Definition: deprecated_function_map.h:75
FEValues< dim >
Subscriptor
Definition: subscriptor.h:62
Timer
Definition: timer.h:119
internal::TriangulationImplementation::n_active_cells
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12625
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
PreconditionSSOR::initialize
void initialize(const MatrixType &A, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
QTrapez
Definition: quadrature_lib.h:126
Physics::Elasticity::Kinematics::w
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition: dof_tools_sparsity.cc:63
VectorizedArray::exp
VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5368
VectorTools::interpolate_boundary_values
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
Algorithms::Events::initial
const Event initial
Definition: event.cc:65
Timer::start
void start()
Definition: timer.cc:177
Tensor< 1, dim >
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
SolutionTransfer
Definition: solution_transfer.h:340
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
KellyErrorEstimator::estimate
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandlerType &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)
DynamicSparsityPattern
Definition: dynamic_sparsity_pattern.h:323
SparsityPattern
Definition: sparsity_pattern.h:865
std::abs
inline ::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5450
Function::value
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, AffineConstraints< number > &constraints)
Definition: dof_tools_constraints.cc:1725
GridGenerator::merge_triangulations
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)
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:66
QGauss
Definition: quadrature_lib.h:40
LocalIntegrators::L2::L2
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:171
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandlerType &)
LogStream::depth_console
unsigned int depth_console(const unsigned int n)
Definition: logstream.cc:349
value
static const bool value
Definition: dof_tools_constraints.cc:433
GridGenerator::hyper_ball
void hyper_ball(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
AffineConstraints< double >
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
DataOutBase::eps
@ eps
Definition: data_out_base.h:1582
std::sqrt
inline ::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5412
VectorTools::H1_seminorm
@ H1_seminorm
Definition: vector_tools_common.h:165
VectorizedArray::sqrt
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5412
GridGenerator::hyper_cube
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
std::pow
inline ::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &x, const Number p)
Definition: vectorization.h:5428
Utilities::MPI::min
T min(const T &t, const MPI_Comm &mpi_communicator)
Point< dim >
Functions::ZeroFunction
Definition: function.h:513
ParameterHandler
Definition: parameter_handler.h:845
PreconditionSSOR
Definition: precondition.h:665
FullMatrix< double >
DataPostprocessor::get_needed_update_flags
virtual UpdateFlags get_needed_update_flags() const =0
Function
Definition: function.h:151
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
std::log
inline ::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5390
SolverControl
Definition: solver_control.h:67
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
DataPostprocessor::get_names
virtual std::vector< std::string > get_names() const =0
first
Point< 2 > first
Definition: grid_out.cc:4352
DataOut< dim >
DoFHandler::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:935
Vector< double >
GridRefinement::refine
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
Definition: grid_refinement.cc:41
ConvergenceTable
Definition: convergence_table.h:63
Patterns::Double
Definition: patterns.h:293
DataPostprocessor::get_data_component_interpretation
virtual std::vector< DataComponentInterpretation::DataComponentInterpretation > get_data_component_interpretation() const
Definition: data_postprocessor.cc:48
DataOut_DoFData::add_data_vector
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation=std::vector< DataComponentInterpretation::DataComponentInterpretation >())
Definition: data_out_dof_data.h:1090
Patterns::Integer
Definition: patterns.h:190