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