clamped_shell_with_arclength_cont.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Driver function for a simple test shell problem:
31 //Calculate the deformation of an elastic tube approximated
32 //using Kirchoff--Love shell theory
33 
34 //Standard system includes
35 #include <iostream>
36 #include <fstream>
37 #include <cmath>
38 #include <typeinfo>
39 #include <algorithm>
40 #include <cstdio>
41 
42 //Include files from the finite-element library
43 #include "generic.h"
44 #include "shell.h"
45 #include "meshes/rectangular_quadmesh.h"
46 
47 using namespace std;
48 
49 using namespace oomph;
50 
51 //========================================================================
52 /// Global variables that represent physical properties
53 //========================================================================
54 namespace Global_Physical_Variables
55 {
56 
57  /// Prescribed position of control point
58  double Prescribed_y = 1.0;
59 
60  /// \short Pointer to pressure load (stored in Data so it can
61  /// become an unknown in the problem when displacement control is used
62  Data* Pext_data_pt;
63 
64  /// \short Return a reference to the external pressure
65  /// load on the elastic tube.
66  double external_pressure()
67  {return (*Pext_data_pt->value_pt(0))*pow(0.05,3)/12.0;}
68 
69 
70  /// Load function, normal pressure loading
71  void press_load(const Vector<double> &xi,
72  const Vector<double> &x,
73  const Vector<double> &N,
74  Vector<double>& load)
75  {
76  for(unsigned i=0;i<3;i++)
77  {
78  load[i] = external_pressure()*N[i];
79  }
80  }
81 
82 }
83 
84 //========================================================================
85 /// A 2D Mesh class. The tube wall is represented by two Lagrangian
86 /// coordinates that correspond to z and theta in cylindrical polars.
87 /// The required mesh is therefore a 2D mesh and is therefore inherited
88 /// from the generic RectangularQuadMesh
89 //=======================================================================
90 template <class ELEMENT>
91 class ShellMesh : public virtual RectangularQuadMesh<ELEMENT>,
92  public virtual SolidMesh
93 {
94 public:
95 
96  ///Constructor for the mesh
97  ShellMesh(const unsigned &nx, const unsigned &ny,
98  const double &lx, const double &ly);
99 
100  /// \short In all elastic problems, the nodes must be assigned an undeformed,
101  /// or reference, position, corresponding to the stress-free state
102  /// of the elastic body. This function assigns the undeformed position
103  /// for the nodes on the elastic tube
104  void assign_undeformed_positions(GeomObject* const &undeformed_midplane_pt);
105 
106 };
107 
108 
109 
110 
111 
112 //=======================================================================
113 /// Mesh constructor
114 /// Argument list:
115 /// nx : number of elements in the axial direction
116 /// ny : number of elements in the azimuthal direction
117 /// lx : length in the axial direction
118 /// ly : length in theta direction
119 //=======================================================================
120 template <class ELEMENT>
121 ShellMesh<ELEMENT>::ShellMesh(const unsigned &nx,
122  const unsigned &ny,
123  const double &lx,
124  const double &ly) :
125  RectangularQuadMesh<ELEMENT>(nx,ny,lx,ly)
126 {
127  //Find out how many nodes there are
128  unsigned n_node = nnode();
129 
130  //Now in this case it is the Lagrangian coordinates that we want to set,
131  //so we have to loop over all nodes and set them to the Eulerian
132  //coordinates that are set by the generic mesh generator
133  for(unsigned i=0;i<n_node;i++)
134  {
135  node_pt(i)->xi(0) = node_pt(i)->x(0);
136  node_pt(i)->xi(1) = node_pt(i)->x(1);
137  }
138 
139 
140  //Assign gradients, etc for the Lagrangian coordinates of
141  //hermite-type elements
142 
143  //Read out number of position dofs
144  unsigned n_position_type = finite_element_pt(0)->nnodal_position_type();
145 
146  //If this is greater than 1 set the slopes, which are the distances between
147  //nodes. If the spacing were non-uniform, this part would be more difficult
148  if(n_position_type > 1)
149  {
150  double xstep = (this->Xmax - this->Xmin)/((this->Np-1)*this->Nx);
151  double ystep = (this->Ymax - this->Ymin)/((this->Np-1)*this->Ny);
152  for(unsigned n=0;n<n_node;n++)
153  {
154  //The factor 0.5 is because our reference element has length 2.0
155  node_pt(n)->xi_gen(1,0) = 0.5*xstep;
156  node_pt(n)->xi_gen(2,1) = 0.5*ystep;
157  }
158  }
159 }
160 
161 
162 //=======================================================================
163 /// Set the undeformed coordinates of the nodes
164 //=======================================================================
165 template <class ELEMENT>
167  GeomObject* const &undeformed_midplane_pt)
168 {
169  //Find out how many nodes there are
170  unsigned n_node = nnode();
171 
172  //Loop over all the nodes
173  for(unsigned n=0;n<n_node;n++)
174  {
175  //Get the Lagrangian coordinates
176  Vector<double> xi(2);
177  xi[0] = node_pt(n)->xi(0);
178  xi[1] = node_pt(n)->xi(1);
179 
180  //Assign memory for values of derivatives, etc
181  Vector<double> R(3);
182  DenseMatrix<double> a(2,3);
183  RankThreeTensor<double> dadxi(2,2,3);
184 
185  //Get the geometrical information from the geometric object
186  undeformed_midplane_pt->d2position(xi,R,a,dadxi);
187 
188  //Loop over coordinate directions
189  for(unsigned i=0;i<3;i++)
190  {
191  //Set the position
192  node_pt(n)->x_gen(0,i) = R[i];
193 
194  //Set the derivative wrt Lagrangian coordinates
195  //Note that we need to scale by the length of each element here!!
196  node_pt(n)->x_gen(1,i) = 0.5*a(0,i)*((this->Xmax - this->Xmin)/this->Nx);
197  node_pt(n)->x_gen(2,i) = 0.5*a(1,i)*((this->Ymax - this->Ymin)/this->Ny);
198 
199  //Set the mixed derivative
200  //(symmetric so doesn't matter which one we use)
201  node_pt(n)->x_gen(3,i) = 0.25*dadxi(0,1,i);
202  }
203  }
204 }
205 
206 
207 
208 // //========================================================================
209 // /// A 2D mesh class for the problem.
210 // /// The tube is represented by two Lagrangian coordinates that correspond
211 // /// to z and theta in cylindrical polars. The required mesh is therefore a
212 // /// 2D mesh and is therefore inherited from the generic RectangularQuadMesh
213 // //=======================================================================
214 // template <class ELEMENT>
215 // class ShellMesh : public virtual RectangularQuadMesh<ELEMENT>,
216 // public virtual SolidMesh
217 // {
218 // public:
219 // //Constructor for the mesh
220 // ShellMesh(const unsigned &nx, const unsigned &ny,
221 // const double &lx, const double &ly);
222 
223 // //In all elastic problems, the nodes must be assigned an undeformed,
224 // //or reference, position, corresponding to the stress-free state
225 // //of the elastic body. This function assigns the undeformed position
226 // //for the nodes on the elastic tube
227 // void assign_undeformed_positions(GeomObject* const &undeformed_midplane_pt);
228 // };
229 
230 // //Define the mesh constructor
231 // //Argument list:
232 // // nx : number of elements in the axial direction
233 // // ny : number of elements in the azimuthal direction
234 // // lx : length in the axial direction
235 // // ly : length in theta direction
236 // template <class ELEMENT>
237 // ShellMesh<ELEMENT>::ShellMesh(const unsigned &nx,
238 // const unsigned &ny,
239 // const double &lx, const double &ly) :
240 // RectangularQuadMesh<ELEMENT>(nx,ny,lx,ly)
241 // {
242 // //Find out how many nodes there are
243 // unsigned Nnode = nnode();
244 
245 // //Now in this case it is the Lagrangian coordinates that we want to set,
246 // //so we have to loop over all nodes and set them to the Eulerian
247 // //coordinates that are set by the generic mesh generator
248 // for(unsigned i=0;i<Nnode;i++)
249 // {
250 // node_pt(i)->xi(0) = node_pt(i)->x(0);
251 // node_pt(i)->xi(1) = node_pt(i)->x(1);
252 // }
253 
254 // //Assign gradients, etc for the Lagrangian coordinates of
255 // //hermite-type elements
256 
257 // //Read out number of position dofs
258 // unsigned Nposition_type = finite_element_pt(0)->nnodal_position_type();
259 // std::cout << Nposition_type << std::endl;
260 // //If this is greater than 1 set the slopes, which are the distances between
261 // //nodes. If the spacing were non-uniform, this part would be more difficult
262 // if(Nposition_type > 1)
263 // {
264 // double xstep = (this->Xmax - this->Xmin)/((this->Np-1)*this->Nx);
265 // double ystep = (this->Ymax - this->Ymin)/((this->Np-1)*this->Ny);
266 // for(unsigned n=0;n<Nnode;n++)
267 // {
268 // //The factor 0.5 is because our reference element has length 2.0
269 // node_pt(n)->xi_gen(1,0) = 0.5*xstep;
270 // node_pt(n)->xi_gen(2,1) = 0.5*ystep;
271 // }
272 // }
273 // }
274 
275 // //Set the undeformed values of the nodes
276 // template <class ELEMENT>
277 // void ShellMesh<ELEMENT>::
278 // assign_undeformed_positions(GeomObject* const &undeformed_midplane_pt)
279 // {
280 // //Find out how many nodes there are
281 // unsigned n_node = nnode();
282 // //Loop over all the nodes
283 // for(unsigned n=0;n<n_node;n++)
284 // {
285 // //Get the lagrangian coordinates
286 // Vector<double> xi(2);
287 // xi[0] = node_pt(n)->xi(0); xi[1] = node_pt(n)->xi(1);
288 
289 // //Assign memory for values of derivatives, etc
290 // Vector<double> R(3);
291 // DenseMatrix<double> a(2,3);
292 // RankThreeTensor<double> dadxi(2,2,3);
293 
294 // //Get the geometrical information from the geometric object
295 // undeformed_midplane_pt->d2position(xi,R,a,dadxi);
296 
297 // //Loop over coordinate directions
298 // for(unsigned i=0;i<3;i++)
299 // {
300 // //Set the position
301 // node_pt(n)->x_gen(0,i) = R[i];
302 // //Set the derivative wrt Lagrangian coordinates
303 // //Note that we need to scale by the length of each element here!!
304 // node_pt(n)->x_gen(1,i) = 0.5*a(0,i)*((this->Xmax - this->Xmin)/this->Nx);
305 // node_pt(n)->x_gen(2,i) = 0.5*a(1,i)*((this->Ymax - this->Ymin)/this->Ny);
306 // //Set the mixed derivative
307 // //(symmetric so doesn't matter which one we use)
308 // node_pt(n)->x_gen(3,i) = 0.25*dadxi(0,1,i);
309 // }
310 // }
311 // }
312 
313 
314 
315 //======================================================================
316 //Problem class to solve the deformation of an elastic tube
317 //=====================================================================
318 template<class ELEMENT>
319 class ShellProblem : public Problem
320 {
321 
322 public:
323 
324  /// Constructor
325  ShellProblem(const unsigned &nx, const unsigned &ny,
326  const double &lx, const double &ly);
327 
328  /// Destructor: delete mesh, geometric object
330  {
331  delete Problem::mesh_pt();
332  delete Undeformed_midplane_pt;
333  }
334 
335 
336  /// Overload Access function for the mesh
338  {return dynamic_cast<ShellMesh<ELEMENT>*>(Problem::mesh_pt());}
339 
340  /// Actions after solve empty
342 
343  /// Actions before solve empty
345 
346  //A self_test function
347  void solve();
348 
349 private:
350 
351  /// Pointer to GeomObject that specifies the undeformed midplane
352  GeomObject* Undeformed_midplane_pt;
353 
354  /// First trace node
355  Node* Trace_node_pt;
356 
357  /// Second trace node
358  Node* Trace_node2_pt;
359 
360 };
361 
362 
363 
364 //======================================================================
365 /// Constructor
366 //======================================================================
367 template<class ELEMENT>
368 ShellProblem<ELEMENT>::ShellProblem(const unsigned &nx, const unsigned &ny,
369  const double &lx, const double &ly)
370 {
371  //Suppress the warning about empty mesh_level_timestepper_stuff
372  Mesh::Suppress_warning_about_empty_mesh_level_time_stepper_function=true;
373 
374  //Use the continuation timestepper
375  Use_continuation_timestepper = true;
376 
377  //Create the undeformed midplane object
378  Undeformed_midplane_pt = new EllipticalTube(1.0,1.0);
379 
380  //Now create the mesh
381  Problem::mesh_pt() = new ShellMesh<ELEMENT>(nx,ny,lx,ly);
382 
383  //Set the undeformed positions in the mesh
384  mesh_pt()->assign_undeformed_positions(Undeformed_midplane_pt);
385 
386  //Reorder the elements, since I know what's best for them....
387  mesh_pt()->element_reorder();
388 
389  //Apply boundary conditions to the ends of the tube
390  unsigned n_ends = mesh_pt()->nboundary_node(1);
391  //Loop over the node
392  for(unsigned i=0;i<n_ends;i++)
393  {
394  //Pin in the axial direction (prevents rigid body motions)
395  mesh_pt()->boundary_node_pt(1,i)->pin_position(2);
396  mesh_pt()->boundary_node_pt(3,i)->pin_position(2);
397  //Derived conditions
398  mesh_pt()->boundary_node_pt(1,i)->pin_position(2,2);
399  mesh_pt()->boundary_node_pt(3,i)->pin_position(2,2);
400 
401  //------------------CLAMPING CONDITIONS----------------------
402  //------Pin positions in the transverse directions-----------
403  // Comment these out to get the ring case
404  mesh_pt()->boundary_node_pt(1,i)->pin_position(0);
405  mesh_pt()->boundary_node_pt(3,i)->pin_position(0);
406  //Derived conditions
407  mesh_pt()->boundary_node_pt(1,i)->pin_position(2,0);
408  mesh_pt()->boundary_node_pt(3,i)->pin_position(2,0);
409 
410  mesh_pt()->boundary_node_pt(1,i)->pin_position(1);
411  mesh_pt()->boundary_node_pt(3,i)->pin_position(1);
412  //Derived conditions
413  mesh_pt()->boundary_node_pt(1,i)->pin_position(2,1);
414  mesh_pt()->boundary_node_pt(3,i)->pin_position(2,1);
415  //----------------------------------------------------------
416 
417  // Set the axial gradients of the transverse coordinates to be
418  // zero --- need to be enforced for ring or tube buckling
419  //Pin dx/dz and dy/dz
420  mesh_pt()->boundary_node_pt(1,i)->pin_position(1,0);
421  mesh_pt()->boundary_node_pt(1,i)->pin_position(1,1);
422  mesh_pt()->boundary_node_pt(3,i)->pin_position(1,0);
423  mesh_pt()->boundary_node_pt(3,i)->pin_position(1,1);
424  //Derived conditions
425  mesh_pt()->boundary_node_pt(1,i)->pin_position(3,0);
426  mesh_pt()->boundary_node_pt(1,i)->pin_position(3,1);
427  mesh_pt()->boundary_node_pt(3,i)->pin_position(3,0);
428  mesh_pt()->boundary_node_pt(3,i)->pin_position(3,1);
429  }
430 
431  //Now loop over the sides and apply symmetry conditions
432  unsigned n_side = mesh_pt()->nboundary_node(0);
433  for(unsigned i=0;i<n_side;i++)
434  {
435  //At the side where theta is 0, pin in the y direction
436  mesh_pt()->boundary_node_pt(0,i)->pin_position(1);
437  //Derived condition
438  mesh_pt()->boundary_node_pt(0,i)->pin_position(1,1);
439  //Pin dx/dtheta and dz/dtheta
440  mesh_pt()->boundary_node_pt(0,i)->pin_position(2,0);
441  mesh_pt()->boundary_node_pt(0,i)->pin_position(2,2);
442  //Pin the mixed derivative
443  mesh_pt()->boundary_node_pt(0,i)->pin_position(3,0);
444  mesh_pt()->boundary_node_pt(0,i)->pin_position(3,2);
445 
446  //At the side when theta is 0.5pi pin in the x direction
447  mesh_pt()->boundary_node_pt(2,i)->pin_position(0);
448  //Derived condition
449  mesh_pt()->boundary_node_pt(2,i)->pin_position(1,0);
450  //Pin dy/dtheta and dz/dtheta
451  mesh_pt()->boundary_node_pt(2,i)->pin_position(2,1);
452  mesh_pt()->boundary_node_pt(2,i)->pin_position(2,2);
453  //Pin the mixed derivative
454  mesh_pt()->boundary_node_pt(2,i)->pin_position(3,1);
455  mesh_pt()->boundary_node_pt(2,i)->pin_position(3,2);
456 
457  //Set an initial kick to make sure that we hop onto the
458  //non-axisymmetric branch
459  if((i>1) && (i<n_side-1))
460  {
461  mesh_pt()->boundary_node_pt(0,i)->x(0) += 0.05;
462  mesh_pt()->boundary_node_pt(2,i)->x(1) -= 0.1;
463  }
464  }
465 
466 
467  // Setup displacement control
468  //---------------------------
469 
470 
471 
472 // //Setup displacement control
473 // //Fix the displacement at the mid-point of the tube in the "vertical"
474 // //(y) direction.
475 // //Set the displacement control element (located halfway along the tube)
476 // Disp_ctl_element_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(3*Ny-1));
477 // //The midpoint of the tube is located exactly half-way along the element
478 // Vector<double> s(2); s[0] = 1.0; s[1] = 0.0; //s[1] = 0.5
479 // //Fix the displacement at this point in the y (1) direction
480 // Disp_ctl_element_pt->fix_displacement_for_displacement_control(s,1);
481 // //Set the pointer to the prescribed position
482 // Disp_ctl_element_pt->prescribed_position_pt() = &Prescribed_y;
483 
484 
485 
486  // Choose element in which displacement control is applied: This
487  // one is located halfway along the tube)
488  SolidFiniteElement* controlled_element_pt=
489  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(3*ny-1));
490 
491  // Fix the displacement in the y (1) direction...
492  unsigned controlled_direction=1;
493 
494  // Local coordinate of the control point within the controlled element
495  Vector<double> s_displ_control(2);
496  s_displ_control[0]=1.0;
497  s_displ_control[1]=0.0;
498 
499  // Pointer to displacement control element
500  DisplacementControlElement* displ_control_el_pt;
501 
502  // Build displacement control element
503  displ_control_el_pt=
504  new DisplacementControlElement(controlled_element_pt,
505  s_displ_control,
506  controlled_direction,
508 
509  // The constructor of the DisplacementControlElement has created
510  // a new Data object whose one-and-only value contains the
511  // adjustable load: Use this Data object in the load function:
512  Global_Physical_Variables::Pext_data_pt=displ_control_el_pt->
513  displacement_control_load_pt();
514 
515  // Add the displacement-control element to the mesh
516  mesh_pt()->add_element_pt(displ_control_el_pt);
517 
518 
519 
520  // Complete build of shell elements
521  //---------------------------------
522 
523  //Find number of shell elements in mesh
524  unsigned n_element = nx*ny;
525 
526  //Explicit pointer to first element in the mesh
527  ELEMENT* first_el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(0));
528 
529  //Loop over the elements
530  for(unsigned e=0;e<n_element;e++)
531  {
532  //Cast to a shell element
533  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
534 
535  //Set the load function
536  el_pt->load_vector_fct_pt() = & Global_Physical_Variables::press_load;
537 
538  //Set the undeformed surface
539  el_pt->undeformed_midplane_pt() = Undeformed_midplane_pt;
540 
541  //The external pressure is external data for all elements
542  el_pt->add_external_data(Global_Physical_Variables::Pext_data_pt);
543 
544  //Pre-compute the second derivatives wrt Lagrangian coordinates
545  //for the first element only
546  if(e==0)
547  {
548  el_pt->pre_compute_d2shape_lagrangian_at_knots();
549  }
550 
551  //Otherwise set the values to be the same as those in the first element
552  //this is OK because the Lagrangian mesh is uniform.
553  else
554  {
555  el_pt->set_dshape_lagrangian_stored_from_element(first_el_pt);
556  }
557  }
558 
559  //Set pointers to two trace nodes, used for output
560  Trace_node_pt = mesh_pt()->finite_element_pt(2*ny-1)->node_pt(3);
561  Trace_node2_pt = mesh_pt()->finite_element_pt(ny)->node_pt(1);
562 
563  // Do equation numbering
564  cout << std::endl;
565  cout << "------------------DISPLACEMENT CONTROL--------------------"
566  << std::endl;
567  cout << "# of dofs " << assign_eqn_numbers() << std::endl;
568  cout << "----------------------------------------------------------"
569  << std::endl;
570  cout << std::endl;
571 
572 }
573 
574 
575 //================================================================
576 // /Define the solve function, disp ctl and then continuation
577 //================================================================
578 template<class ELEMENT>
580 {
581 
582  //Increase the maximum number of Newton iterations.
583  //Finding the first buckled solution requires a large(ish) number
584  //of Newton steps -- shells are just a bit twitchy
585  Max_newton_iterations = 40;
586 
587  //Open an output trace file
588  ofstream trace("trace.dat");
589 
590  //Change in displacemenet
591  double disp_incr = 0.05;
592 
593  //Gradually compress the tube by decreasing the value of the prescribed
594  //position from 1 to zero in steps of 0.05 initially and then 0.1
595  for(unsigned i=1;i<13;i++)
596  {
597  //By the time we reach the second time round increase the incremenet
598  if(i==3) {disp_incr = 0.1;}
599  //Reduce prescribed y by our chosen increment
601 
602  cout << std::endl << "Increasing displacement: Prescribed_y is "
604 
605  // Solve
606  newton_solve();
607 
608  //Output the pressure
609  trace << Global_Physical_Variables::external_pressure()/(pow(0.05,3)/12.0)
610  << " "
611  //Position of first trace node
612  << Trace_node_pt->x(0) << " " << Trace_node_pt->x(1) << " "
613  //Position of second trace node
614  << Trace_node2_pt->x(0) << " " << Trace_node2_pt->x(1) << std::endl;
615  }
616 
617  //Close the trace file
618  trace.close();
619 
620  //Output the tube shape in the most strongly collapsed configuration
621  ofstream file("final_shape.dat");
622  mesh_pt()->output(file,5);
623  file.close();
624 
625 
626  //Switch from displacement control to arc-length continuation and
627  //trace back up the solution branch
628 
629  //Now pin the external pressure
631 
632  //Re-assign the equation numbers
633  cout << std::endl;
634  cout << "-----------------ARC-LENGTH CONTINUATION --------------"
635  << std::endl;
636  cout << "# of dofs " << assign_eqn_numbers() << std::endl;
637  cout << "-------------------------------------------------------"
638  << std::endl;
639  cout << std::endl;
640 
641  //Set the maximum number of Newton iterations to something more reasonable
642  Max_newton_iterations=6;
643 
644  //Set the desired number of Newton iterations per arc-length step
645  Desired_newton_iterations_ds=3;
646 
647  //Set the proportion of the arc length
648  Desired_proportion_of_arc_length = 0.2;
649 
650  //Check for any bifurcations using the sign of the determinant of the Jacobian
651  Bifurcation_detection = true;
652 
653  //Set an initial value for the step size
654  double ds = -0.5;
655 
656  //Open a different trace file
657  trace.open("trace_disp.dat");
658  //Take fifteen continuation steps
659  for(unsigned i=0;i<15;i++)
660  {
661  ds = arc_length_step_solve(
663 
664  //Output the pressure
665  trace << Global_Physical_Variables::external_pressure()/(pow(0.05,3)/12.0)
666  << " "
667  //Position of first trace node
668  << Trace_node_pt->x(0) << " " << Trace_node_pt->x(1) << " "
669  //Position of second trace node
670  << Trace_node2_pt->x(0) << " " << Trace_node2_pt->x(1) << std::endl;
671  }
672 
673  //Close the trace file
674  trace.close();
675 
676 }
677 
678 
679 //====================================================================
680 /// Driver
681 //====================================================================
682 int main()
683 {
684 
685  //Length of domain
686  double L = 10.0;
687  double L_phi=0.5*MathematicalConstants::Pi;
688 
689  //Set up the problem
691  problem(5,5,L,L_phi);
692 
693  //Solve the problem
694  problem.solve();
695 }
696 
697 
698 
699 
700 
701 
void actions_before_newton_solve()
Actions before solve empty.
void assign_undeformed_positions(GeomObject *const &undeformed_midplane_pt)
In all elastic problems, the nodes must be assigned an undeformed, or reference, position, corresponding to the stress-free state of the elastic body. This function assigns the undeformed position for the nodes on the elastic tube.
Data * Pext_data_pt
Pointer to pressure load (stored in Data so it can become an unknown in the problem when displacement...
void press_load(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Load function, normal pressure loading.
Node * Trace_node_pt
First trace node.
GeomObject * Undeformed_midplane_pt
Pointer to GeomObject that specifies the undeformed midplane.
void actions_after_newton_solve()
Actions after solve empty.
double external_pressure()
Return a reference to the external pressure load on the elastic tube.
~ShellProblem()
Destructor: delete mesh, geometric object.
double Prescribed_y
Prescribed position of control point.
int main()
Driver.
ShellProblem(const unsigned &nx, const unsigned &ny, const double &lx, const double &ly)
Constructor.
ShellMesh(const unsigned &nx, const unsigned &ny, const double &lx, const double &ly)
Constructor for the mesh.
ShellMesh< ELEMENT > * mesh_pt()
Overload Access function for the mesh.
Node * Trace_node2_pt
Second trace node.