clamped_shell.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  /// Perturbation pressure
65  double Pcos=1.0;
66 
67 
68  /// \short Return a reference to the external pressure
69  /// load on the elastic tube.
71  {return (*Pext_data_pt->value_pt(0))*pow(0.05,3)/12.0;}
72 
73 
74  /// Load function, normal pressure loading
75  void press_load(const Vector<double> &xi,
76  const Vector<double> &x,
77  const Vector<double> &N,
78  Vector<double>& load)
79  {
80  //std::cout << N[0] << " " << N[1] << " " << N[2] << std::endl;
81  //std::cout << xi[0] << " " << xi[1] << std::endl;
82  for(unsigned i=0;i<3;i++)
83  {
84  load[i] = (external_pressure() -
85  Pcos*pow(0.05,3)/12.0*cos(2.0*xi[1]))*N[i];
86  }
87  }
88 
89 }
90 
91 //========================================================================
92 /// A 2D Mesh class. The tube wall is represented by two Lagrangian
93 /// coordinates that correspond to z and theta in cylindrical polars.
94 /// The required mesh is therefore a 2D mesh and is therefore inherited
95 /// from the generic RectangularQuadMesh
96 //=======================================================================
97 template <class ELEMENT>
98 class ShellMesh : public virtual RectangularQuadMesh<ELEMENT>,
99  public virtual SolidMesh
100 {
101 public:
102 
103  ///Constructor for the mesh
104  ShellMesh(const unsigned &nx, const unsigned &ny,
105  const double &lx, const double &ly);
106 
107  /// \short In all elastic problems, the nodes must be assigned an undeformed,
108  /// or reference, position, corresponding to the stress-free state
109  /// of the elastic body. This function assigns the undeformed position
110  /// for the nodes on the elastic tube
111  void assign_undeformed_positions(GeomObject* const &undeformed_midplane_pt);
112 
113 };
114 
115 
116 
117 
118 
119 //=======================================================================
120 /// Mesh constructor
121 /// Argument list:
122 /// nx : number of elements in the axial direction
123 /// ny : number of elements in the azimuthal direction
124 /// lx : length in the axial direction
125 /// ly : length in theta direction
126 //=======================================================================
127 template <class ELEMENT>
128 ShellMesh<ELEMENT>::ShellMesh(const unsigned &nx,
129  const unsigned &ny,
130  const double &lx,
131  const double &ly) :
132  RectangularQuadMesh<ELEMENT>(nx,ny,lx,ly)
133 {
134  //Find out how many nodes there are
135  unsigned n_node = nnode();
136 
137  //Now in this case it is the Lagrangian coordinates that we want to set,
138  //so we have to loop over all nodes and set them to the Eulerian
139  //coordinates that are set by the generic mesh generator
140  for(unsigned i=0;i<n_node;i++)
141  {
142  node_pt(i)->xi(0) = node_pt(i)->x(0);
143  node_pt(i)->xi(1) = node_pt(i)->x(1);
144  }
145 
146 
147  //Assign gradients, etc for the Lagrangian coordinates of
148  //hermite-type elements
149 
150  //Read out number of position dofs
151  unsigned n_position_type = finite_element_pt(0)->nnodal_position_type();
152 
153  //If this is greater than 1 set the slopes, which are the distances between
154  //nodes. If the spacing were non-uniform, this part would be more difficult
155  if(n_position_type > 1)
156  {
157  double xstep = (this->Xmax - this->Xmin)/((this->Np-1)*this->Nx);
158  double ystep = (this->Ymax - this->Ymin)/((this->Np-1)*this->Ny);
159  for(unsigned n=0;n<n_node;n++)
160  {
161  //The factor 0.5 is because our reference element has length 2.0
162  node_pt(n)->xi_gen(1,0) = 0.5*xstep;
163  node_pt(n)->xi_gen(2,1) = 0.5*ystep;
164  }
165  }
166 }
167 
168 
169 //=======================================================================
170 /// Set the undeformed coordinates of the nodes
171 //=======================================================================
172 template <class ELEMENT>
174  GeomObject* const &undeformed_midplane_pt)
175 {
176  //Find out how many nodes there are
177  unsigned n_node = nnode();
178 
179  //Loop over all the nodes
180  for(unsigned n=0;n<n_node;n++)
181  {
182  //Get the Lagrangian coordinates
183  Vector<double> xi(2);
184  xi[0] = node_pt(n)->xi(0);
185  xi[1] = node_pt(n)->xi(1);
186 
187  //Assign memory for values of derivatives, etc
188  Vector<double> R(3);
189  DenseMatrix<double> a(2,3);
190  RankThreeTensor<double> dadxi(2,2,3);
191 
192  //Get the geometrical information from the geometric object
193  undeformed_midplane_pt->d2position(xi,R,a,dadxi);
194 
195  //Loop over coordinate directions
196  for(unsigned i=0;i<3;i++)
197  {
198  //Set the position
199  node_pt(n)->x_gen(0,i) = R[i];
200 
201  //Set the derivative wrt Lagrangian coordinates
202  //Note that we need to scale by the length of each element here!!
203  node_pt(n)->x_gen(1,i) = 0.5*a(0,i)*((this->Xmax - this->Xmin)/this->Nx);
204  node_pt(n)->x_gen(2,i) = 0.5*a(1,i)*((this->Ymax - this->Ymin)/this->Ny);
205 
206  //Set the mixed derivative
207  //(symmetric so doesn't matter which one we use)
208  node_pt(n)->x_gen(3,i) = 0.25*dadxi(0,1,i);
209  }
210  }
211 }
212 
213 
214 //======================================================================
215 //Problem class to solve the deformation of an elastic tube
216 //=====================================================================
217 template<class ELEMENT>
218 class ShellProblem : public Problem
219 {
220 
221 public:
222 
223  /// Constructor
224  ShellProblem(const unsigned &nx, const unsigned &ny,
225  const double &lx, const double &ly);
226 
227  /// Overload Access function for the mesh
229  {return dynamic_cast<ShellMesh<ELEMENT>*>(Problem::mesh_pt());}
230 
231  /// Actions after solve empty
233 
234  /// Actions before solve empty
236 
237  //A self_test function
238  void solve();
239 
240 private:
241 
242  /// Pointer to GeomObject that specifies the undeformed midplane
244 
245  /// First trace node
247 
248  /// Second trace node
250 
251 };
252 
253 
254 
255 //======================================================================
256 /// Constructor
257 //======================================================================
258 template<class ELEMENT>
259 ShellProblem<ELEMENT>::ShellProblem(const unsigned &nx, const unsigned &ny,
260  const double &lx, const double &ly)
261 {
262  //Create the undeformed midplane object
263  Undeformed_midplane_pt = new EllipticalTube(1.0,1.0);
264 
265  //Now create the mesh
266  Problem::mesh_pt() = new ShellMesh<ELEMENT>(nx,ny,lx,ly);
267 
268  //Set the undeformed positions in the mesh
269  mesh_pt()->assign_undeformed_positions(Undeformed_midplane_pt);
270 
271  //Reorder the elements, since I know what's best for them....
272  mesh_pt()->element_reorder();
273 
274  //Apply boundary conditions to the ends of the tube
275  unsigned n_ends = mesh_pt()->nboundary_node(1);
276  //Loop over the node
277  for(unsigned i=0;i<n_ends;i++)
278  {
279  //Pin in the axial direction (prevents rigid body motions)
280  mesh_pt()->boundary_node_pt(1,i)->pin_position(2);
281  mesh_pt()->boundary_node_pt(3,i)->pin_position(2);
282  //Derived conditions
283  mesh_pt()->boundary_node_pt(1,i)->pin_position(2,2);
284  mesh_pt()->boundary_node_pt(3,i)->pin_position(2,2);
285 
286  //------------------CLAMPING CONDITIONS----------------------
287  //------Pin positions in the transverse directions-----------
288  // Comment these out to get the ring case
289  mesh_pt()->boundary_node_pt(1,i)->pin_position(0);
290  mesh_pt()->boundary_node_pt(3,i)->pin_position(0);
291  //Derived conditions
292  mesh_pt()->boundary_node_pt(1,i)->pin_position(2,0);
293  mesh_pt()->boundary_node_pt(3,i)->pin_position(2,0);
294 
295  mesh_pt()->boundary_node_pt(1,i)->pin_position(1);
296  mesh_pt()->boundary_node_pt(3,i)->pin_position(1);
297  //Derived conditions
298  mesh_pt()->boundary_node_pt(1,i)->pin_position(2,1);
299  mesh_pt()->boundary_node_pt(3,i)->pin_position(2,1);
300  //----------------------------------------------------------
301 
302  // Set the axial gradients of the transverse coordinates to be
303  // zero --- need to be enforced for ring or tube buckling
304  //Pin dx/dz and dy/dz
305  mesh_pt()->boundary_node_pt(1,i)->pin_position(1,0);
306  mesh_pt()->boundary_node_pt(1,i)->pin_position(1,1);
307  mesh_pt()->boundary_node_pt(3,i)->pin_position(1,0);
308  mesh_pt()->boundary_node_pt(3,i)->pin_position(1,1);
309  //Derived conditions
310  mesh_pt()->boundary_node_pt(1,i)->pin_position(3,0);
311  mesh_pt()->boundary_node_pt(1,i)->pin_position(3,1);
312  mesh_pt()->boundary_node_pt(3,i)->pin_position(3,0);
313  mesh_pt()->boundary_node_pt(3,i)->pin_position(3,1);
314  }
315 
316  //Now loop over the sides and apply symmetry conditions
317  unsigned n_side = mesh_pt()->nboundary_node(0);
318  for(unsigned i=0;i<n_side;i++)
319  {
320  //At the side where theta is 0, pin in the y direction
321  mesh_pt()->boundary_node_pt(0,i)->pin_position(1);
322  //Derived condition
323  mesh_pt()->boundary_node_pt(0,i)->pin_position(1,1);
324  //Pin dx/dtheta and dz/dtheta
325  mesh_pt()->boundary_node_pt(0,i)->pin_position(2,0);
326  mesh_pt()->boundary_node_pt(0,i)->pin_position(2,2);
327  //Pin the mixed derivative
328  mesh_pt()->boundary_node_pt(0,i)->pin_position(3,0);
329  mesh_pt()->boundary_node_pt(0,i)->pin_position(3,2);
330 
331  //At the side when theta is 0.5pi pin in the x direction
332  mesh_pt()->boundary_node_pt(2,i)->pin_position(0);
333  //Derived condition
334  mesh_pt()->boundary_node_pt(2,i)->pin_position(1,0);
335  //Pin dy/dtheta and dz/dtheta
336  mesh_pt()->boundary_node_pt(2,i)->pin_position(2,1);
337  mesh_pt()->boundary_node_pt(2,i)->pin_position(2,2);
338  //Pin the mixed derivative
339  mesh_pt()->boundary_node_pt(2,i)->pin_position(3,1);
340  mesh_pt()->boundary_node_pt(2,i)->pin_position(3,2);
341 
342 // //Set an initial kick to make sure that we hop onto the
343 // //non-axisymmetric branch
344 // if((i>1) && (i<n_side-1))
345 // {
346 // mesh_pt()->boundary_node_pt(0,i)->x(0) += 0.05;
347 // mesh_pt()->boundary_node_pt(2,i)->x(1) -= 0.1;
348 // }
349  }
350 
351 
352  // Setup displacement control
353  //---------------------------
354 
355 
356 
357 // //Setup displacement control
358 // //Fix the displacement at the mid-point of the tube in the "vertical"
359 // //(y) direction.
360 // //Set the displacement control element (located halfway along the tube)
361 // Disp_ctl_element_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(3*Ny-1));
362 // //The midpoint of the tube is located exactly half-way along the element
363 // Vector<double> s(2); s[0] = 1.0; s[1] = 0.0; //s[1] = 0.5
364 // //Fix the displacement at this point in the y (1) direction
365 // Disp_ctl_element_pt->fix_displacement_for_displacement_control(s,1);
366 // //Set the pointer to the prescribed position
367 // Disp_ctl_element_pt->prescribed_position_pt() = &Prescribed_y;
368 
369 
370 
371  // Choose element in which displacement control is applied: This
372  // one is located about halfway along the tube -- remember that
373  // we've renumbered the elements!
374  unsigned nel_ctrl=0;
375  Vector<double> s_displ_control(2);
376 
377  // Even/odd number of elements in axial direction
378  if (nx%2==1)
379  {
380  nel_ctrl=unsigned(floor(0.5*double(nx))+1.0)*ny-1;
381  s_displ_control[0]=0.0;
382  s_displ_control[1]=1.0;
383  }
384  else
385  {
386  nel_ctrl=unsigned(floor(0.5*double(nx))+1.0)*ny-1;
387  s_displ_control[0]=-1.0;
388  s_displ_control[1]=1.0;
389  }
390 
391  // Controlled element
392  SolidFiniteElement* controlled_element_pt=
393  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(nel_ctrl));
394 
395 
396  // Fix the displacement in the y (1) direction...
397  unsigned controlled_direction=1;
398 
399  // Pointer to displacement control element
400  DisplacementControlElement* displ_control_el_pt;
401 
402  // Build displacement control element
403  displ_control_el_pt=
404  new DisplacementControlElement(controlled_element_pt,
405  s_displ_control,
406  controlled_direction,
408 
409 
410  // Doc control point
411  Vector<double> xi(2);
412  Vector<double> x(3);
413  controlled_element_pt->interpolated_xi(s_displ_control,xi);
414  controlled_element_pt->interpolated_x(s_displ_control,x);
415  std::cout << std::endl;
416  std::cout << "Controlled element: " << nel_ctrl << std::endl;
417  std::cout << "Displacement control applied at xi = ("
418  << xi[0] << ", " << xi[1] << ")" << std::endl;
419  std::cout << "Corresponding to x = ("
420  << x[0] << ", " << x[1] << ", " << x[2] << ")" << std::endl;
421 
422 
423  // The constructor of the DisplacementControlElement has created
424  // a new Data object whose one-and-only value contains the
425  // adjustable load: Use this Data object in the load function:
426  Global_Physical_Variables::Pext_data_pt=displ_control_el_pt->
427  displacement_control_load_pt();
428 
429  // Add the displacement-control element to the mesh
430  mesh_pt()->add_element_pt(displ_control_el_pt);
431 
432 
433 
434  // Complete build of shell elements
435  //---------------------------------
436 
437  //Find number of shell elements in mesh
438  unsigned n_element = nx*ny;
439 
440  //Explicit pointer to first element in the mesh
441  ELEMENT* first_el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(0));
442 
443  //Loop over the elements
444  for(unsigned e=0;e<n_element;e++)
445  {
446  //Cast to a shell element
447  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
448 
449  //Set the load function
450  el_pt->load_vector_fct_pt() = & Global_Physical_Variables::press_load;
451 
452  //Set the undeformed surface
453  el_pt->undeformed_midplane_pt() = Undeformed_midplane_pt;
454 
455  //The external pressure is external data for all elements
456  el_pt->add_external_data(Global_Physical_Variables::Pext_data_pt);
457 
458  //Pre-compute the second derivatives wrt Lagrangian coordinates
459  //for the first element only
460  if(e==0)
461  {
462  el_pt->pre_compute_d2shape_lagrangian_at_knots();
463  }
464 
465  //Otherwise set the values to be the same as those in the first element
466  //this is OK because the Lagrangian mesh is uniform.
467  else
468  {
469  el_pt->set_dshape_lagrangian_stored_from_element(first_el_pt);
470  }
471  }
472 
473  //Set pointers to two trace nodes, used for output
474  Trace_node_pt = mesh_pt()->finite_element_pt(2*ny-1)->node_pt(3);
475  Trace_node2_pt = mesh_pt()->finite_element_pt(ny)->node_pt(1);
476 
477  // Do equation numbering
478  cout << std::endl;
479  cout << "# of dofs " << assign_eqn_numbers() << std::endl;
480  cout << std::endl;
481 
482 }
483 
484 
485 //================================================================
486 // /Define the solve function, disp ctl and then continuation
487 //================================================================
488 template<class ELEMENT>
490 {
491 
492  //Increase the maximum number of Newton iterations.
493  //Finding the first buckled solution requires a large(ish) number
494  //of Newton steps -- shells are just a bit twitchy
495  Max_newton_iterations = 40;
496  Max_residuals=1.0e6;
497 
498 
499  //Open an output trace file
500  ofstream trace("trace.dat");
501 
502 
503  //Gradually compress the tube by decreasing the value of the prescribed
504  //position
505  for(unsigned i=1;i<11;i++)
506  {
507 
509 
510  cout << std::endl << "Increasing displacement: Prescribed_y is "
512 
513  // Solve
514  newton_solve();
515 
516 
517  //Output the pressure (on the bending scale)
518  trace << Global_Physical_Variables::external_pressure()/(pow(0.05,3)/12.0)
519  << " "
520  //Position of first trace node
521  << Trace_node_pt->x(0) << " " << Trace_node_pt->x(1) << " "
522  //Position of second trace node
523  << Trace_node2_pt->x(0) << " " << Trace_node2_pt->x(1) << std::endl;
524 
525  // Reset perturbation
527  }
528 
529  //Close the trace file
530  trace.close();
531 
532  //Output the tube shape in the most strongly collapsed configuration
533  ofstream file("final_shape.dat");
534  mesh_pt()->output(file,5);
535  file.close();
536 
537 
538 }
539 
540 
541 //====================================================================
542 /// Driver
543 //====================================================================
544 int main()
545 {
546 
547  //Length of domain
548  double L = 10.0;
549  double L_phi=0.5*MathematicalConstants::Pi;
550 
551  //Set up the problem
553  problem(5,3,L,L_phi);
554 
555  //Solve the problem
556  problem.solve();
557 }
558 
559 
560 
561 
562 
563 
void actions_before_newton_solve()
Actions before solve empty.
double Pcos
Perturbation pressure.
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.
int main()
Driver.
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.
double Prescribed_y
Prescribed position of control point.
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.