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