Reference documentation for deal.II version 9.5.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
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 override;
241 *   virtual Tensor<1,dim> gradient (const Point<dim> &pto, const unsigned int component = 0) const override;
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 override;
296 *  
297 *   virtual
298 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
299 *   get_data_component_interpretation () const override;
300 *  
301 *   virtual UpdateFlags get_needed_update_flags () const override;
302 *   };
303 *  
304 *  
305 *   template <int dim>
306 *   ComputeMultiplier<dim>::ComputeMultiplier (double pe): p(pe)
307 *   {}
308 *  
309 *  
310 *   template <int dim>
311 *   void ComputeMultiplier<dim>::compute_derived_quantities_scalar(
312 *   const std::vector< double > & /*uh*/,
313 *   const std::vector< Tensor< 1, dim > > &duh,
314 *   const std::vector< Tensor< 2, dim > > & /*dduh*/,
315 *   const std::vector< Point< dim > > & /* normals*/,
316 *   const std::vector< Point< dim > > & /*evaluation_points*/,
317 *   std::vector< Vector< double > > &computed_quantities ) const
318 *   {
319 *   const unsigned int n_quadrature_points = duh.size();
320 *  
321 *   for (unsigned int q=0; q<n_quadrature_points; ++q)
322 *   {
323 *   long double sqrGrad=duh[q]* duh[q]; //squared norm of the gradient
324 *   long double exponent=(p-2.0)/2*std::log(sqrGrad);
325 *   computed_quantities[q](0) = std::sqrt(sqrGrad); // norm of the gradient
326 *   computed_quantities[q](1)= std::exp(exponent); // multiplier
327 *   }
328 *   }
329 *  
330 *  
331 *  
332 *  
333 *  
334 *   template <int dim>
335 *   std::vector<std::string>
336 *   ComputeMultiplier<dim>::get_names() const
337 *   {
338 *   std::vector<std::string> solution_names;
339 *   solution_names.push_back ("Gradient norm");
340 *   solution_names.push_back ("Lagrange multiplier");
341 *   return solution_names;
342 *   }
343 *  
344 *  
345 *   template <int dim>
346 *   UpdateFlags
347 *   ComputeMultiplier<dim>::get_needed_update_flags () const
348 *   {
349 *   return update_gradients;
350 *   }
351 *  
352 *  
353 *  
354 *   template <int dim>
355 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
356 *   ComputeMultiplier<dim>:: get_data_component_interpretation () const
357 *   {
358 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
359 *   interpretation;
360 * @endcode
361 *
362 * norm of the gradient
363 *
364 * @code
365 *   interpretation.push_back (DataComponentInterpretation::component_is_scalar);
366 * @endcode
367 *
368 * Lagrange multiplier
369 *
370 * @code
371 *   interpretation.push_back (DataComponentInterpretation::component_is_scalar);
372 *   return interpretation;
373 *   }
374 *  
375 *  
376 *  
377 *  
378 *  
379 * @endcode
380 *
381 * ***************************************************************************************
382 *
383 * @code
384 *   template <int dim>
385 *   class ElastoplasticTorsion
386 *   {
387 *   public:
388 *   ElastoplasticTorsion (ParameterHandler &);
389 *   ~ElastoplasticTorsion ();
390 *   void run ();
391 *  
392 *   private:
393 *   void setup_system (const bool initial_step);
394 *   void assemble_system ();
395 *   bool solve (const int inner_it);
396 *   void init_mesh ();
397 *   void refine_mesh ();
398 *   void set_boundary_values ();
399 *   double phi (const double alpha) const;
400 *   bool checkWolfe(double &alpha, double &phi_alpha) const;
401 *   bool determine_step_length (const int inner_it);
402 *   void print_it_message (const int counter, bool ks);
403 *   void output_results (unsigned int refinement) const;
404 *   void format_convergence_tables();
405 *   void process_solution (const unsigned int cycle);
406 *   void process_multiplier (const unsigned int cycle,const int iter,double time);
407 *   double dual_error () const;
408 *   double dual_infty_error () const;
409 *   double W (double Du2) const;
410 *   double Wp (double Du2) const;
411 *   double G (double Du2) const;
412 *  
413 *  
414 *  
415 *   ParameterHandler &prm;
417 *   DoFHandler<dim> dof_handler;
418 *   AffineConstraints<double> hanging_node_constraints;
419 *   SparsityPattern sparsity_pattern;
420 *   SparseMatrix<double> system_matrix;
421 *   ConvergenceTable convergence_table;
422 *   ConvergenceTable dual_convergence_table;
423 *   Vector<double> present_solution;
424 *   Vector<double> newton_update;
425 *   Vector<double> system_rhs;
426 *   Vector<double> grad_norm;
428 *  
429 *  
430 *   double step_length,phi_zero,phi_alpha,phip,phip_zero;
431 *   double old_step,old_phi_zero,old_phip;
432 *   double L2_error;
433 *   double H1_error;
434 *   double Linfty_error;
435 *   double dual_L1_error;
436 *   double dual_L_infty_error;
437 *   FE_Q<dim> fe;
438 *   double p;
439 *   double line_search_tolerence; // c_1 in Nocedal & Wright
440 *   unsigned int dir_id;
441 *   std::string elements;
442 *   std::string Method;
443 *  
444 *   };
445 *  
446 *   /*******************************************************************************************/
447 * @endcode
448 *
449 * Boundary condition
450 *
451
452 *
453 *
454 * @code
455 *   template <int dim>
456 *   class BoundaryValues : public Function<dim>
457 *   {
458 *   public:
459 *   BoundaryValues () : Function<dim>() {}
460 *  
461 *   virtual double value (const Point<dim> &p,
462 *   const unsigned int component = 0) const override;
463 *   };
464 *  
465 *  
466 *   template <int dim>
467 *   double BoundaryValues<dim>::value (const Point<dim> &/*pto*/,
468 *   const unsigned int /*component*/) const
469 *   {
470 * @endcode
471 *
472 * could be anything else (theory works provided |Dg|_infty < 1/2)
473 *
474 * @code
475 *   return 0.0;
476 *  
477 *   /* A challenging BC leading to overdetermined problems
478 *   it is regulated by the parameter 0<eta<1.
479 *   eta closer to 1 leads to more difficult problems.
480 *  
481 *   double pii=numbers::PI;
482 *   double theta=std::atan2(p[1],p[0])+pii;
483 *   double eta=0.9;
484 *  
485 *   if (theta <= 0.5)
486 *   return eta*(theta*theta);
487 *   else if ((theta >0.5) & (theta<= pii-0.5))
488 *   return eta*(theta-0.25);
489 *   else if ((theta>pii-0.5)&(theta<= pii+0.5))
490 *   return eta*(pii-0.75-(theta-(pii-0.5))*(theta-(pii+0.5)));
491 *   else if ((theta>pii+0.5)&(theta<= 2*pii-0.5))
492 *   return eta*((2*pii-theta)-0.25);
493 *   else
494 *   return eta*((theta-2*pii)*(theta-2*pii) );*/
495 *   }
496 *  
497 *  
498 *  
499 *   /******************************************************************************/
500 * @endcode
501 *
502 * Right-Hand Side
503 *
504 * @code
505 *   template <int dim>
506 *   class RightHandSide : public Function<dim>
507 *   {
508 *   public:
509 *   RightHandSide () : Function<dim>() {}
510 *   virtual double value (const Point<dim> &p,
511 *   const unsigned int component = 0) const override;
512 *   };
513 *  
514 *   template <int dim>
515 *   double RightHandSide<dim>::value (const Point<dim> &/*p*/,
516 *   const unsigned int /*component*/) const
517 *   {
518 * @endcode
519 *
520 * set to constant = 4, for which explicit solution to compare exists
521 * could be anything
522 *
523 * @code
524 *   double return_value = 4.0;
525 *   return return_value;
526 *   }
527 *  
528 *  
529 *  
530 *   /*******************************************************************/
531 * @endcode
532 *
533 * The ElastoplasticTorsion class implementation
534 *
535
536 *
537 * Constructor of the class
538 *
539 * @code
540 *   template <int dim>
541 *   ElastoplasticTorsion<dim>::ElastoplasticTorsion (ParameterHandler &param):
542 *   prm(param),
543 *   dof_handler (triangulation),
544 *   L2_error(1.0),
545 *   H1_error(1.0),
546 *   Linfty_error(1.0),
547 *   dual_L1_error(1.0),
548 *   dual_L_infty_error(1.0),
549 *   fe(2)
550 *   {
551 *   prm.enter_subsection ("Global Parameters");
552 *   p=prm.get_double("p");
553 *   prm.leave_subsection ();
554 *   prm.enter_subsection ("Algorithm Parameters");
555 *   line_search_tolerence=prm.get_double("line_search_tolerence");
556 *   dir_id=prm.get_integer("Descent_direction");
557 *   prm.leave_subsection ();
558 *   if (fe.degree==1)
559 *   elements="P1";
560 *   else elements="P2";
561 *  
562 *   if (dir_id==0)
563 *   Method="Precond";
564 *   else
565 *   Method="Newton";
566 *   }
567 *  
568 *  
569 *  
570 *   template <int dim>
571 *   ElastoplasticTorsion<dim>::~ElastoplasticTorsion ()
572 *   {
573 *   dof_handler.clear ();
574 *   }
575 *  
576 *   /*****************************************************************************************/
577 * @endcode
578 *
579 * print iteration message
580 *
581
582 *
583 *
584 * @code
585 *   template <int dim>
586 *   void ElastoplasticTorsion<dim>::print_it_message (const int counter, bool ks)
587 *   {
588 *   if (ks)
589 *   {
590 *   process_solution (counter);
591 *   std::cout << "iteration="<< counter+1 << " J(u_h)= "<< phi_zero << ", H1 error: "
592 *   << H1_error <<", W0-1,infty error: "<< Linfty_error<< " J'(u_h)(w)= "<< phip
593 *   << ", |J'(u_h)|= "<< system_rhs.l2_norm()<<std::endl;
594 *   }
595 *   else
596 *   {
597 *   std::cout << "iteration= " << counter+1 << " J(u_h)= "
598 *   << phi_alpha << " J'(u_h)= "<< phip<<std::endl;
599 *   }
600 *   }
601 *  
602 *  
603 *   /*****************************************************************************************/
604 * @endcode
605 *
606 * Convergence Tables
607 *
608
609 *
610 *
611
612 *
613 *
614 * @code
615 *   /*************************************************************/
616 * @endcode
617 *
618 * formating
619 *
620
621 *
622 *
623 * @code
624 *   template <int dim>
625 *   void ElastoplasticTorsion<dim>::format_convergence_tables()
626 *   {
627 *   convergence_table.set_precision("L2", 3);
628 *   convergence_table.set_precision("H1", 3);
629 *   convergence_table.set_precision("Linfty", 3);
630 *   convergence_table.set_precision("function value", 3);
631 *   convergence_table.set_precision("derivative", 3);
632 *   dual_convergence_table.set_precision("dual_L1", 3);
633 *   dual_convergence_table.set_precision("dual_Linfty", 3);
634 *   dual_convergence_table.set_precision("L2", 3);
635 *   dual_convergence_table.set_precision("H1", 3);
636 *   dual_convergence_table.set_precision("Linfty", 3);
637 *   convergence_table.set_scientific("L2", true);
638 *   convergence_table.set_scientific("H1", true);
639 *   convergence_table.set_scientific("Linfty", true);
640 *   convergence_table.set_scientific("function value", true);
641 *   convergence_table.set_scientific("derivative", true);
642 *   dual_convergence_table.set_scientific("dual_L1", true);
643 *   dual_convergence_table.set_scientific("dual_Linfty", true);
644 *   dual_convergence_table.set_scientific("L2", true);
645 *   dual_convergence_table.set_scientific("H1", true);
646 *   dual_convergence_table.set_scientific("Linfty", true);
647 *  
648 *   }
649 *  
650 *   /****************************************/
651 * @endcode
652 *
653 * fill-in entry for the solution
654 *
655 * @code
656 *   template <int dim>
657 *   void ElastoplasticTorsion<dim>::process_solution (const unsigned int it)
658 *   {
659 *   Vector<float> difference_per_cell (triangulation.n_active_cells());
660 *  
661 * @endcode
662 *
663 * compute L2 error (save to difference_per_cell)
664 *
665 * @code
666 *   VectorTools::integrate_difference (dof_handler,present_solution,
667 *   Solution<dim>(),difference_per_cell,QGauss<dim>(3),VectorTools::L2_norm);
668 *   L2_error = difference_per_cell.l2_norm();
669 *  
670 * @endcode
671 *
672 * compute H1 error (save to difference_per_cell)
673 *
674 * @code
675 *   VectorTools::integrate_difference (dof_handler,present_solution,Solution<dim>(),
676 *   difference_per_cell,QGauss<dim>(3),VectorTools::H1_seminorm);
677 *   H1_error = difference_per_cell.l2_norm();
678 *  
679 * @endcode
680 *
681 * compute W1infty error (save to difference_per_cell)
682 *
683 * @code
684 *   const QTrapezoid<1> q_trapez;
685 *   const QIterated<dim> q_iterated (q_trapez, 5);
686 *   VectorTools::integrate_difference (dof_handler,present_solution,Solution<dim>(),
687 *   difference_per_cell,q_iterated,VectorTools::W1infty_seminorm);
688 *   Linfty_error = difference_per_cell.linfty_norm();
689 *  
690 *  
691 *   convergence_table.add_value("cycle", it);
692 *   convergence_table.add_value("p", p);
693 *   convergence_table.add_value("L2", L2_error);
694 *   convergence_table.add_value("H1", H1_error);
695 *   convergence_table.add_value("Linfty", Linfty_error);
696 *   convergence_table.add_value("function value", phi_alpha);
697 *   convergence_table.add_value("derivative", phip);
698 *   }
699 *  
700 *  
701 *   /***************************************/
702 * @endcode
703 *
704 * fill-in entry for the multiplier
705 *
706 * @code
707 *   template <int dim>
708 *   void ElastoplasticTorsion<dim>::process_multiplier (const unsigned int cycle, const int iter,double time)
709 *   {
710 *   const unsigned int n_active_cells=triangulation.n_active_cells();
711 *   const unsigned int n_dofs=dof_handler.n_dofs();
712 *   dual_L1_error=dual_error();
713 *   dual_L_infty_error=dual_infty_error();
714 *  
715 *  
716 *   dual_convergence_table.add_value("cycle", cycle);
717 *   dual_convergence_table.add_value("p", p);
718 *   dual_convergence_table.add_value("iteration_number", iter);
719 *   dual_convergence_table.add_value("cpu_time", time);
720 *   dual_convergence_table.add_value("cells", n_active_cells);
721 *   dual_convergence_table.add_value("dofs", n_dofs);
722 *   dual_convergence_table.add_value("L2", L2_error);
723 *   dual_convergence_table.add_value("H1", H1_error);
724 *   dual_convergence_table.add_value("Linfty", Linfty_error);
725 *   dual_convergence_table.add_value("dual_L1", dual_L1_error);
726 *   dual_convergence_table.add_value("dual_Linfty", dual_L_infty_error);
727 *  
728 *   }
729 *  
730 *  
731 *  
732 *  
733 *   /****************************************************************************************/
734 * @endcode
735 *
736 * ElastoplasticTorsion::setup_system
737 * unchanged from @ref step_15 "step-15"
738 *
739
740 *
741 *
742 * @code
743 *   template <int dim>
744 *   void ElastoplasticTorsion<dim>::setup_system (const bool initial_step)
745 *   {
746 *   if (initial_step)
747 *   {
748 *   dof_handler.distribute_dofs (fe);
749 *   present_solution.reinit (dof_handler.n_dofs());
750 *   grad_norm.reinit (dof_handler.n_dofs());
751 *   lambda.reinit (dof_handler.n_dofs());
752 *  
753 *   hanging_node_constraints.clear ();
755 *   hanging_node_constraints);
756 *   hanging_node_constraints.close ();
757 *   }
758 *  
759 *  
760 * @endcode
761 *
762 * The remaining parts of the function
763 *
764
765 *
766 *
767 * @code
768 *   newton_update.reinit (dof_handler.n_dofs());
769 *   system_rhs.reinit (dof_handler.n_dofs());
770 *   DynamicSparsityPattern c_sparsity(dof_handler.n_dofs());
771 *   DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
772 *   hanging_node_constraints.condense (c_sparsity);
773 *   sparsity_pattern.copy_from(c_sparsity);
774 *   system_matrix.reinit (sparsity_pattern);
775 *   }
776 *  
777 *   /***************************************************************************************/
778 *   /* the coeffcients W, W' and G defining the problem.
779 *  
780 *   Min_u \int W(|Du|^2) dx
781 *  
782 *   They must be consistent as G(s)=W'(s)+2s W''(s) for any s>0.
783 *   recall that they receive the SQUARED gradient. */
784 *  
785 *   template <int dim>
786 *   double ElastoplasticTorsion<dim>::W (double Du2) const
787 *   {
788 *   return Du2;
789 *   }
790 *  
791 *   template <int dim>
792 *   double ElastoplasticTorsion<dim>::Wp (double /*Du2*/) const
793 *   {
794 *   return 1.0;
795 *   }
796 *  
797 *   template <int dim>
798 *   double ElastoplasticTorsion<dim>::G (double /*Du2*/) const
799 *   {
800 *   return 1.0;
801 *   }
802 *   /***************************************************************************************/
803 *  
804 *   template <int dim>
805 *   void ElastoplasticTorsion<dim>::assemble_system ()
806 *   {
807 *  
808 *  
809 *   const QGauss<dim> quadrature_formula(3);
810 *   const RightHandSide<dim> right_hand_side;
811 *   system_matrix = 0;
812 *   system_rhs = 0;
813 *  
814 *   FEValues<dim> fe_values (fe, quadrature_formula,
815 *   update_gradients |
816 *   update_values |
818 *   update_JxW_values);
819 *  
820 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
821 *   const unsigned int n_q_points = quadrature_formula.size();
822 *  
823 *   FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
824 *   Vector<double> cell_rhs (dofs_per_cell);
825 *  
826 *   std::vector<Tensor<1, dim> > old_solution_gradients(n_q_points);
827 *   std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
828 *  
829 *  
831 *   cell = dof_handler.begin_active(),
832 *   endc = dof_handler.end();
833 *   for (; cell!=endc; ++cell)
834 *   {
835 *   cell_matrix = 0;
836 *   cell_rhs = 0;
837 *  
838 *   fe_values.reinit (cell);
839 *   fe_values.get_function_gradients(present_solution,
840 *   old_solution_gradients);
841 *  
842 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
843 *   {
844 *   long double coeff=0.0;
845 *   long double a=old_solution_gradients[q_point] * old_solution_gradients[q_point];
846 *   long double exponent=(p-2.0)/2*std::log(a);
847 *   coeff= std::exp( exponent);
848 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
849 *   {
850 *   for (unsigned int j=0; j<dofs_per_cell; ++j)
851 *   {
852 *   if (dir_id==1)
853 *   {
854 *   cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * fe_values.shape_grad(j, q_point)
855 *   * (G(a)+(p-1.0)*coeff) * fe_values.JxW(q_point);
856 *   }
857 *   else
858 *   {
859 *   cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * fe_values.shape_grad(j, q_point)
860 *   * (Wp(a)+coeff)
861 *   * fe_values.JxW(q_point);
862 *   }
863 *   }
864 *  
865 *   cell_rhs(i) -= ( fe_values.shape_grad(i, q_point)
866 *   * old_solution_gradients[q_point]
867 *   * (Wp(a)+coeff)
868 *   -right_hand_side.value(fe_values.quadrature_point(q_point))
869 *   *fe_values.shape_value(i, q_point)
870 *   )
871 *   * fe_values.JxW(q_point);
872 *   }
873 *   }
874 *  
875 *   cell->get_dof_indices (local_dof_indices);
876 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
877 *   {
878 *   for (unsigned int j=0; j<dofs_per_cell; ++j)
879 *   system_matrix.add (local_dof_indices[i],
880 *   local_dof_indices[j],
881 *   cell_matrix(i,j));
882 *  
883 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
884 *   }
885 *   }
886 *  
887 *   hanging_node_constraints.condense (system_matrix);
888 *   hanging_node_constraints.condense (system_rhs);
889 *  
890 *   std::map<types::global_dof_index,double> boundary_values;
892 *   0,
894 *   boundary_values);
895 *   MatrixTools::apply_boundary_values (boundary_values,
896 *   system_matrix,
897 *   newton_update,
898 *   system_rhs);
899 *   }
900 *  
901 *  
902 *  
903 *  
904 *   /********************************** Refine Mesh ****************************************/
905 * @endcode
906 *
907 * unchanged from @ref step_15 "step-15"
908 *
909
910 *
911 *
912 * @code
913 *   template <int dim>
914 *   void ElastoplasticTorsion<dim>::refine_mesh ()
915 *   {
916 *   using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
917 *  
918 *   Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
919 *   KellyErrorEstimator<dim>::estimate (dof_handler,
920 *   QGauss<dim-1>(3),
921 *   FunctionMap(),
922 *   present_solution,
923 *   estimated_error_per_cell);
924 *  
925 *   prm.enter_subsection ("Mesh & Refinement Parameters");
926 *   const double top_fraction=prm.get_double("top_fraction_of_cells");
927 *   const double bottom_fraction=prm.get_double("bottom_fraction_of_cells");
928 *   prm.leave_subsection ();
930 *   estimated_error_per_cell,
931 *   top_fraction, bottom_fraction);
932 *  
933 *   triangulation.prepare_coarsening_and_refinement ();
934 *   SolutionTransfer<dim> solution_transfer(dof_handler);
935 *   solution_transfer.prepare_for_coarsening_and_refinement(present_solution);
936 *   triangulation.execute_coarsening_and_refinement();
937 *   dof_handler.distribute_dofs(fe);
938 *   Vector<double> tmp(dof_handler.n_dofs());
939 *   solution_transfer.interpolate(present_solution, tmp);
940 *   present_solution = tmp;
941 *   set_boundary_values ();
942 *   hanging_node_constraints.clear();
943 *  
945 *   hanging_node_constraints);
946 *   hanging_node_constraints.close();
947 *   hanging_node_constraints.distribute (present_solution);
948 *   setup_system (false);
949 *   }
950 *  
951 *  
952 *   /*******************************************************************************************/
953 * @endcode
954 *
955 * Dump the norm of the gradient and the lagrange multiplier in vtu format for visualization
956 *
957 * @code
958 *   template <int dim>
959 *   void ElastoplasticTorsion<dim>::output_results (unsigned int counter) const
960 *   {
961 * @endcode
962 *
963 * multiplier object contains both |Du| and lambda.
964 *
965 * @code
966 *   ComputeMultiplier<dim> multiplier(p);
967 *   DataOut<dim> data_out;
968 *  
969 *   data_out.attach_dof_handler (dof_handler);
970 *   data_out.add_data_vector (present_solution, "solution");
971 *   data_out.add_data_vector (present_solution, multiplier);
972 *   data_out.build_patches ();
973 *   std::ostringstream p_str;
974 *   p_str << p<<"-cycle-"<<counter;
975 *   std::string str = p_str.str();
976 *   const std::string filename = "solution-" + str+".vtu";
977 *   std::ofstream output (filename.c_str());
978 *   data_out.write_vtu (output);
979 *   }
980 *  
981 *   /********************************************************************************************/
982 * @endcode
983 *
984 * unchanged from @ref step_15 "step-15"
985 *
986 * @code
987 *   template <int dim>
988 *   void ElastoplasticTorsion<dim>::set_boundary_values ()
989 *   {
990 *   std::map<types::global_dof_index, double> boundary_values;
992 *   0,
993 *   BoundaryValues<dim>(),
994 *   boundary_values);
995 *   for (std::map<types::global_dof_index, double>::const_iterator
996 *   bp = boundary_values.begin();
997 *   bp != boundary_values.end(); ++bp)
998 *   present_solution(bp->first) = bp->second;
999 *   }
1000 *  
1001 *  
1002 *   /****************************************************************************************/
1003 * @endcode
1004 *
1005 * COMPUTE @f$\phi(\alpha)=J_p(u_h+\alpha w)@f$
1006 *
1007 * @code
1008 *   template <int dim>
1009 *   double ElastoplasticTorsion<dim>::phi (const double alpha) const
1010 *   {
1011 *   double obj = 0.0;
1012 *   const RightHandSide<dim> right_hand_side;
1013 *   Vector<double> evaluation_point (dof_handler.n_dofs());
1014 *   evaluation_point = present_solution; // copy of u_h
1015 *   evaluation_point.add (alpha, newton_update); // u_{n+1}=u_n+alpha w_n
1016 *  
1017 *   const QGauss<dim> quadrature_formula(3);
1018 *   FEValues<dim> fe_values (fe, quadrature_formula,
1019 *   update_gradients |
1020 *   update_values |
1022 *   update_JxW_values);
1023 *  
1024 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
1025 *   const unsigned int n_q_points = quadrature_formula.size();
1026 *  
1027 *   Vector<double> cell_residual (dofs_per_cell);
1028 *   std::vector<Tensor<1, dim> > gradients(n_q_points);
1029 *   std::vector<double> values(n_q_points);
1030 *  
1031 *  
1032 *   std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1033 *  
1035 *   cell = dof_handler.begin_active(),
1036 *   endc = dof_handler.end();
1037 *   for (; cell!=endc; ++cell)
1038 *   {
1039 *   cell_residual = 0;
1040 *   fe_values.reinit (cell);
1041 *   fe_values.get_function_gradients (evaluation_point, gradients);
1042 *   fe_values.get_function_values (evaluation_point, values);
1043 *  
1044 *  
1045 *   for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
1046 *   {
1047 *   double Du2=gradients[q_point] * gradients[q_point]; // Du2=|Du|^2
1048 *   double penalty;
1049 *   if (Du2<1.0e-10)
1050 *   penalty=0.0;
1051 *   else
1052 *   penalty=std::pow(Du2,p/2.0); // penalty=|Du|^p
1053 *  
1054 * @endcode
1055 *
1056 * obj+= 1/2 W(|Du|^2)+1/p |Du|^p -fu (see (1))
1057 *
1058 * @code
1059 *   obj+=(
1060 *   (0.5*W(Du2)+penalty/p)- right_hand_side.value(fe_values.quadrature_point(q_point))*values[q_point]
1061 *   ) * fe_values.JxW(q_point);
1062 *   }
1063 *  
1064 *   }
1065 *  
1066 *   return obj;
1067 *   }
1068 *  
1069 *  
1070 *   /***************************************************************************************************/
1071 * @endcode
1072 *
1073 * Compute L^1 error norm of Lagrange Multiplier
1074 * with respect to exact solution (cf. Alvarez & Flores, 2015)
1075 *
1076
1077 *
1078 *
1079 * @code
1080 *   template <int dim>
1081 *   double ElastoplasticTorsion<dim>::dual_error () const
1082 *   {
1083 *   double obj = 0.0;
1084 *  
1085 *   const QGauss<dim> quadrature_formula(3);
1086 *   FEValues<dim> fe_values (fe, quadrature_formula,
1087 *   update_gradients |
1089 *   update_JxW_values);
1090 *  
1091 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
1092 *   const unsigned int n_q_points = quadrature_formula.size();
1093 *  
1094 *   Vector<double> cell_residual (dofs_per_cell);
1095 *   std::vector<Tensor<1, dim> > gradients(n_q_points);
1096 *  
1097 *   std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1098 *  
1100 *   cell = dof_handler.begin_active(),
1101 *   endc = dof_handler.end();
1102 *   for (; cell!=endc; ++cell)
1103 *   {
1104 *   cell_residual = 0;
1105 *   fe_values.reinit (cell);
1106 *   fe_values.get_function_gradients (present_solution, gradients);
1107 *  
1108 *   for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
1109 *   {
1110 *   double coeff=gradients[q_point] * gradients[q_point] ;
1111 *   if (coeff<1.0e-15)
1112 *   coeff=0.0;
1113 *   else
1114 *   coeff=std::pow(coeff,(p-2.0)/2.0); // |Du_p|^(p-2)
1115 *  
1116 *   double r=std::sqrt(fe_values.quadrature_point(q_point).square());
1117 *   double exact=0;
1118 *   if (r>0.5)
1119 *   exact= 2*r-1;
1120 *  
1121 *   obj+=( std::abs(coeff-exact) ) * fe_values.JxW(q_point);
1122 *   }
1123 *  
1124 *   }
1125 *  
1126 *   return obj;
1127 *   }
1128 *  
1129 *   /*******************************************************************************************/
1130 * @endcode
1131 *
1132 * Compute L^infinity error norm of Lagrange Multiplier
1133 * with respect to exact solution (cf. Alvarez & Flores, 2015)
1134 *
1135
1136 *
1137 *
1138 * @code
1139 *   template <int dim>
1140 *   double ElastoplasticTorsion<dim>::dual_infty_error () const
1141 *   {
1142 *   double obj = 0.0;
1143 *   const QTrapezoid<1> q_trapez;
1144 *   const QIterated<dim> quadrature_formula (q_trapez, 10);
1145 *  
1146 *   FEValues<dim> fe_values (fe, quadrature_formula,
1147 *   update_gradients |
1149 *  
1150 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
1151 *   const unsigned int n_q_points = quadrature_formula.size();
1152 *  
1153 *   Vector<double> cell_residual (dofs_per_cell);
1154 *   std::vector<Tensor<1, dim> > gradients(n_q_points);
1155 *  
1156 *   std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1157 *  
1159 *   cell = dof_handler.begin_active(),
1160 *   endc = dof_handler.end();
1161 *   for (; cell!=endc; ++cell)
1162 *   {
1163 *   cell_residual = 0;
1164 *   fe_values.reinit (cell);
1165 *   fe_values.get_function_gradients (present_solution, gradients);
1166 *  
1167 *   for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
1168 *   {
1169 *   long double sqdGrad=gradients[q_point] * gradients[q_point] ;
1170 *   double r=std::sqrt(fe_values.quadrature_point(q_point).square());
1171 *   double exact=0;
1172 *   if (r>0.5)
1173 *   exact= 2*r-1.0;
1174 * @endcode
1175 *
1176 * compute |Du|^(p-2) as exp(p-2/2*log(Du^2))
1177 *
1178 * @code
1179 *   long double exponent=(p-2.0)/2*std::log(sqdGrad);
1180 *   long double coeff=std::exp(exponent);
1181 *  
1182 *   if (std::abs(coeff-exact)>obj )
1183 *   obj=std::abs(coeff-exact);
1184 *   }
1185 *  
1186 *   }
1187 *  
1188 *   return obj;
1189 *   }
1190 *  
1191 *   /*****************************************************************************************/
1192 * @endcode
1193 *
1194 * check whether putative step-length satisfies sufficient decrease conditions
1195 *
1196 * @code
1197 *   template <int dim>
1198 *   bool ElastoplasticTorsion<dim>::checkWolfe(double &alpha, double &phi_alpha) const
1199 *   {
1200 *   if (phi_alpha< phi_zero+line_search_tolerence*phip*alpha )
1201 *   return true;
1202 *   else
1203 *   return false;
1204 *   }
1205 *  
1206 *  
1207 *   /*****************************************************************************************/
1208 * @endcode
1209 *
1210 * Find a step-length satisfying sufficient decrease condition by line-search
1211 * uses quadratic interpolation
1212 *
1213
1214 *
1215 *
1216 * @code
1217 *   template <int dim>
1218 *   bool ElastoplasticTorsion<dim>::determine_step_length(const int inner_it)
1219 *   {
1220 *   unsigned int it=0;
1221 *   bool done;
1222 *   double alpha,nalpha;
1223 *   prm.enter_subsection ("Algorithm Parameters");
1224 *   const unsigned int max_LS_it=prm.get_integer("max_LS_it");
1225 *   double init_SL=prm.get_double("init_step_length");
1226 *   prm.leave_subsection ();
1227 *   if (inner_it==0)
1228 *   alpha=init_SL;
1229 *   else
1230 *   {
1231 *   alpha=std::min(1.45*old_step*old_phip/phip,1.0);
1232 *   }
1233 *   phi_alpha=phi(alpha);
1234 *   std::cerr << "Step length=" << alpha << ", Value= " << phi_alpha;
1235 * @endcode
1236 *
1237 * check if step-size satisfies sufficient decrease condition
1238 *
1239 * @code
1240 *   done=checkWolfe(alpha,phi_alpha);
1241 *   if (done)
1242 *   std::cerr << " accepted" << std::endl;
1243 *   else
1244 *   std::cerr << " rejected" ;
1245 *  
1246 *   while ((!done) & (it<max_LS_it))
1247 *   {
1248 * @endcode
1249 *
1250 * new try obtained by quadratic interpolation
1251 *
1252 * @code
1253 *   nalpha=-(phip*alpha*alpha)/(2*(phi_alpha-phi_zero-phip*alpha));
1254 *  
1255 *   if (nalpha<1e-3*alpha || std::abs(nalpha-alpha)/alpha<1e-8)
1256 *   nalpha=alpha/2;
1257 *   else if ( phi_alpha-phi_zero>1e3*std::abs(phi_zero) )
1258 *   nalpha=alpha/10;
1259 *   alpha=nalpha;
1260 *   phi_alpha=phi(alpha);
1261 *   done=checkWolfe(alpha,phi_alpha);
1262 *   if (done)
1263 *   std::cerr << ", finished with steplength= "<< alpha<< ", fcn value= "<< phi_alpha<<std::endl;
1264 *   it=it+1;
1265 *   }
1266 *   if (!done)
1267 *   {
1268 *   std::cerr << ", max. no. of iterations reached wiht steplength= "<< alpha
1269 *   << ", fcn value= "<< phi_alpha<<std::endl;
1270 *   return false;
1271 *   }
1272 *   else
1273 *   {
1274 *   step_length=alpha;
1275 *   return true;
1276 *   }
1277 *  
1278 *   }
1279 *  
1280 *   /**************************************************************************************************/
1281 * @endcode
1282 *
1283 * ElastoplasticTorsion::init_mesh()
1284 *
1285
1286 *
1287 *
1288 * @code
1289 *   template <int dim>
1290 *   void ElastoplasticTorsion<dim>::init_mesh ()
1291 *   {
1292 * @endcode
1293 *
1294 * get parameters
1295 *
1296 * @code
1297 *   prm.enter_subsection ("Mesh & Refinement Parameters");
1298 *   const int domain_id=prm.get_integer("Code for the domain");
1299 *   const int init_ref=prm.get_integer("No of initial refinements");
1300 *   prm.leave_subsection ();
1301 *  
1302 *  
1303 *   if (domain_id==0)
1304 *   {
1305 * @endcode
1306 *
1307 * For the unit disk around the origin
1308 *
1309 * @code
1311 *   static const SphericalManifold<dim> boundary;
1312 *   triangulation.set_manifold (0, boundary);
1313 *   }
1314 *   else if (domain_id==1)
1315 *   {
1316 * @endcode
1317 *
1318 * For the unit square
1319 *
1320 * @code
1322 *   }
1323 *   else if (domain_id==2)
1324 *   {
1325 *   /* For Glowinski's domain
1326 *   ___ ___ __ 1
1327 *   | |__| | __ .8
1328 *   | |
1329 *   | |
1330 *   |__________| __ 0
1331 *  
1332 *   | | | |
1333 *   0 .4 .6 1
1334 *  
1335 *   */
1336 *   Triangulation<dim> tria1;
1337 *   Triangulation<dim> tria2;
1338 *   Triangulation<dim> tria3;
1339 *   Triangulation<dim> tria4;
1340 *   Triangulation<dim> tria5;
1341 *   Triangulation<dim> tria6;
1342 *   GridGenerator::hyper_rectangle(tria1, Point<2>(0.0,0.0), Point<2>(0.4,0.8));
1343 *   GridGenerator::hyper_rectangle(tria2, Point<2>(0.0,0.8), Point<2>(0.4,1.0));
1344 *   GridGenerator::hyper_rectangle(tria3, Point<2>(0.4,0.0), Point<2>(0.6,0.8));
1345 *   GridGenerator::hyper_rectangle(tria4, Point<2>(0.6,0.0), Point<2>(1.0,0.8));
1346 *   GridGenerator::hyper_rectangle(tria5, Point<2>(0.6,0.8), Point<2>(1.0,1.0));
1347 *   GridGenerator::merge_triangulations (tria1, tria2, tria6);
1348 *   GridGenerator::merge_triangulations (tria6, tria3, tria6);
1349 *   GridGenerator::merge_triangulations (tria6, tria4, tria6);
1351 *   }
1352 * @endcode
1353 *
1354 * perform initial refinements
1355 *
1356 * @code
1357 *   triangulation.refine_global(init_ref);
1358 *   }
1359 *  
1360 *   /**************************************************************************************************/
1361 * @endcode
1362 *
1363 * ElastoplasticTorsion::solve(inner_it)
1364 * Performs one inner iteration
1365 *
1366
1367 *
1368 *
1369 * @code
1370 *   template <int dim>
1371 *   bool ElastoplasticTorsion<dim>::solve (const int inner_it)
1372 *   {
1373 *   prm.enter_subsection ("Algorithm Parameters");
1374 *   const unsigned int max_CG_it=prm.get_integer("Max_CG_it");
1375 *   const double CG_tol=prm.get_double("CG_tol");
1376 *   prm.leave_subsection ();
1377 *  
1378 *   SolverControl solver_control (max_CG_it,CG_tol);
1379 *   SolverCG<> solver (solver_control);
1380 *  
1381 *   PreconditionSSOR<> preconditioner;
1382 *   preconditioner.initialize(system_matrix,0.25);
1383 *  
1384 *   solver.solve (system_matrix, newton_update, system_rhs,
1385 *   preconditioner);
1386 *   hanging_node_constraints.distribute (newton_update);
1387 *   /****** save current quantities for line-search **** */
1388 * @endcode
1389 *
1390 * Recall that phi(alpha)=J(u+alpha w)
1391 *
1392 * @code
1393 *   old_step=step_length;
1394 *   old_phi_zero=phi_zero;
1395 *   phi_zero=phi(0); // phi(0)=J(u)
1396 *   old_phip=phip;
1397 *   phip=-1.0*(newton_update*system_rhs); //phi'(0)=J'(u) *w, rhs=-J'(u).
1398 *   if (inner_it==0)
1399 *   phip_zero=phip;
1400 *  
1401 *   if (phip>0) // this should not happen, step back
1402 *   {
1403 *   std::cout << "Not a descent direction!" <<std::endl;
1404 *   present_solution.add (-1.0*step_length, newton_update);
1405 *   step_length=step_length/2;
1406 *   phip=old_phip;
1407 *   return false;
1408 *   }
1409 *   else
1410 *   {
1411 *   if (determine_step_length(inner_it))
1412 *   {
1413 * @endcode
1414 *
1415 * update u_{n+1}=u_n+alpha w_n
1416 *
1417 * @code
1418 *   present_solution.add (step_length, newton_update);
1419 *   return true;
1420 *   }
1421 *   else return false;
1422 *   }
1423 *   }
1424 *  
1425 *  
1426 *  
1427 *   /*************************************************************************************************************/
1428 * @endcode
1429 *
1430 * ElastoplasticTorsion::run
1431 *
1432 * @code
1433 *   template <int dim>
1434 *   void ElastoplasticTorsion<dim>::run ()
1435 *   {
1436 *  
1437 * @endcode
1438 *
1439 * get parameters
1440 *
1441 * @code
1442 *   prm.enter_subsection ("Mesh & Refinement Parameters");
1443 *   const int adapt_ref=prm.get_integer("No of adaptive refinements");
1444 *   prm.leave_subsection ();
1445 *   prm.enter_subsection ("Algorithm Parameters");
1446 *   const int max_inner=prm.get_integer("Max_inner");
1447 *   const double eps=prm.get_double("eps");
1448 *   const double hi_eps=prm.get_double("hi_eps");
1449 *   const int hi_th=prm.get_integer("hi_th");
1450 *   const double init_p=prm.get_double("init_p");
1451 *   const double delta_p=prm.get_double("delta_p");
1452 *   prm.leave_subsection ();
1453 *   prm.enter_subsection ("Global Parameters");
1454 *   bool known_solution=prm.get_bool("known_solution");
1455 *   double actual_p=prm.get_double("p");
1456 *   prm.leave_subsection ();
1457 *   /************************/
1458 *  
1459 * @endcode
1460 *
1461 * init Timer
1462 *
1463 * @code
1464 *   Timer timer;
1465 *   double ptime=0.0;
1466 *   timer.start ();
1467 *  
1468 * @endcode
1469 *
1470 * initalize mesh for the selected domain
1471 *
1472 * @code
1473 *   init_mesh();
1474 *  
1475 * @endcode
1476 *
1477 * setup FE space
1478 *
1479 * @code
1480 *   setup_system (true);
1481 *   set_boundary_values ();
1482 *  
1483 * @endcode
1484 *
1485 * init counters
1486 *
1487 * @code
1488 *   int global_it=0; // Total inner iterations (counting both loops)
1489 *   int cycle=0; // Total outer iterations (counting both loops)
1490 *   int refinement = 0; // Refinements performed (adaptive) = outer iterations 2nd loop
1491 *  
1492 *  
1493 * @endcode
1494 *
1495 * prepare to start first loop
1496 *
1497 * @code
1498 *   p=init_p;
1499 *   bool well_solved=true;
1500 *  
1501 *   /***************************** First loop ***********************************/
1502 *   /****************** Prepare initial condition using increasing p *************************/
1503 *   while (p<actual_p) // outer iteration, increasing p.
1504 *   {
1505 *   std::cout <<"--Preparing initial condition with p="<<p<<" iter.= " << global_it<< " .-- "<< std::endl;
1506 *   timer.restart();
1507 *   for (int inner_iteration=0; inner_iteration<max_inner; ++inner_iteration,++global_it)
1508 *   {
1509 *   assemble_system ();
1510 *   well_solved=solve (inner_iteration);
1511 *   print_it_message (global_it, known_solution);
1512 *   if (
1513 *   ((system_rhs.l2_norm()/std::sqrt(system_rhs.size()) <1e-4) & (cycle<1)) |
1514 *   ((system_rhs.l2_norm()/std::sqrt(system_rhs.size()) <1e-5) & (cycle>=1)) |
1515 *   !well_solved
1516 *   )
1517 *   break;
1518 *   }
1519 *   ptime=timer.cpu_time();
1520 *   if (well_solved)
1521 *   output_results (cycle);
1522 *  
1523 *   if (known_solution)
1524 *   {
1525 *   process_multiplier(cycle,global_it,ptime);
1526 * @endcode
1527 *
1528 * dual_convergence_table.write_tex(dual_error_table_file);
1529 *
1530 * @code
1531 *   }
1532 *   refine_mesh();
1533 *   cycle++;
1534 *   p+=delta_p;
1535 *   }
1536 *   /*************************** first loop finished ********************/
1537 *  
1538 *  
1539 * @endcode
1540 *
1541 * prepare for second loop
1542 *
1543 * @code
1544 *   p=actual_p;
1545 *   well_solved=true;
1546 *  
1547 *  
1548 *   /***************************** Second loop *********************************/
1549 *   /**************************** Solve problem for target p *********************************/
1550 *  
1551 *   std::cout << "============ Solving problem with p=" <<p << " ==================" << std::endl;
1552 *   /***** Outer iteration - refining mesh ******************/
1553 *   while ((cycle<adapt_ref) & well_solved)
1554 *   {
1555 *   timer.restart();
1556 * @endcode
1557 *
1558 * inner iteration
1559 *
1560 * @code
1561 *   for (int inner_iteration=0; inner_iteration<max_inner; ++inner_iteration,++global_it)
1562 *   {
1563 *   assemble_system ();
1564 *   well_solved=solve (inner_iteration);
1565 *   print_it_message (global_it, known_solution);
1566 *  
1567 *   if (
1568 *   ((system_rhs.l2_norm()/std::sqrt(system_rhs.size()) < eps) & (refinement<hi_th)) |
1569 *   (( system_rhs.l2_norm()/ std::sqrt (system_rhs.size()) <hi_eps) | (!well_solved))
1570 *   )
1571 *   break;
1572 *   }
1573 * @endcode
1574 *
1575 * inner iterations finished
1576 *
1577 * @code
1578 *   ptime=timer.cpu_time();
1579 *   if (well_solved)
1580 *   output_results (cycle);
1581 *  
1582 * @endcode
1583 *
1584 * compute and display error, if the explicit solution is known
1585 *
1586 * @code
1587 *   if (known_solution)
1588 *   {
1589 *   process_multiplier(cycle,global_it,ptime);
1590 *   std::cout << "finished with H1 error: " << H1_error << ", dual error (L1): "
1591 *   << dual_L1_error << "dual error (L infty): "<<dual_L_infty_error <<std::endl;
1592 *   }
1593 *  
1594 * @endcode
1595 *
1596 * update counters
1597 *
1598 * @code
1599 *   ++refinement;
1600 *   ++cycle;
1601 * @endcode
1602 *
1603 * refine mesh
1604 *
1605 * @code
1606 *   std::cout << "******** Refined mesh " << cycle << " ********" << std::endl;
1607 *   refine_mesh();
1608 *   }// second loop
1609 *  
1610 * @endcode
1611 *
1612 * write convergence tables to file
1613 *
1614 * @code
1615 *   if (known_solution)
1616 *   {
1617 *   format_convergence_tables();
1618 *   std::string error_filename = "error"+Method+elements+".tex";
1619 *   std::ofstream error_table_file(error_filename.c_str());
1620 *   std::string dual_error_filename = "dual_error"+Method+elements+".tex";
1621 *   std::ofstream dual_error_table_file(dual_error_filename.c_str());
1622 *   convergence_table.write_tex(error_table_file);
1623 *   dual_convergence_table.write_tex(dual_error_table_file);
1624 *   }
1625 *   }//run()
1626 *  
1627 *   }//namespace
1628 *  
1629 *   /**********************************************************************************************/
1630 * @endcode
1631 *
1632 * The main function
1633 *
1634 * @code
1635 *   int main ()
1636 *   {
1637 *   try
1638 *   {
1639 *   using namespace dealii;
1640 *   using namespace nsp;
1641 *   deallog.depth_console (0);
1642 *  
1643 *   ParameterHandler prm;
1644 *   ParameterReader param(prm);
1645 *   param.read_parameters("EPT.prm");
1646 *   ElastoplasticTorsion<2> ElastoplasticTorsionProblem(prm);
1647 *   ElastoplasticTorsionProblem .run ();
1648 *   }
1649 *   catch (std::exception &exc)
1650 *   {
1651 *   std::cerr << std::endl << std::endl
1652 *   << "----------------------------------------------------" << std::endl;
1653 *   std::cerr << "Exception on processing: " << std::endl
1654 *   << exc.what() << std::endl
1655 *   << "Aborting!" << std::endl
1656 *   << "----------------------------------------------------"
1657 *   << std::endl;
1658 *  
1659 *   return 1;
1660 *   }
1661 *   catch (...)
1662 *   {
1663 *   std::cerr << std::endl << std::endl
1664 *   << "----------------------------------------------------"
1665 *   << std::endl;
1666 *   std::cerr << "Unknown exception!" << std::endl
1667 *   << "Aborting!" << std::endl
1668 *   << "----------------------------------------------------"
1669 *   << std::endl;
1670 *   return 1;
1671 *   }
1672 *   return 0;
1673 *   }
1674 * @endcode
1675
1676
1677*/
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:551
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, 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)
unsigned int depth_console(const unsigned int n)
Definition logstream.cc:350
Definition point.h:112
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
Definition timer.h:118
void start()
Definition timer.cc:177
Point< 2 > second
Definition grid_out.cc:4616
Point< 2 > first
Definition grid_out.cc:4615
__global__ void set(Number *val, const Number s, const size_type N)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > 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:439
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=AffineConstraints< number >(), 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.
LogStream deallog
Definition logstream.cc:37
const Event initial
Definition event.cc:65
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:75
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:131
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:472
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition l2.h:160
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 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.)
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=ComponentMask())
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:13833
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
const ::Triangulation< dim, spacedim > & tria