cylinder.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: 1102 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-01-08 08:05:34 +0000 (Fri, 08 Jan 2016) $
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
31 #include <fenv.h>
32 
33 // The oomphlib headers
34 #include "generic.h"
35 #include "axisym_linear_elasticity.h"
36 
37 // The mesh
38 #include "meshes/rectangular_quadmesh.h"
39 
40 using namespace std;
41 
42 using namespace oomph;
43 
44 //===start_of_Global_Parameters_namespace===============================
45 /// Namespace for global parameters
46 //======================================================================
47 namespace Global_Parameters
48 {
49  /// Define Poisson's ratio Nu
50  double Nu = 0.3;
51 
52  /// Define the non-dimensional Young's modulus
53  double E = 1.0;
54 
55  /// Lame parameters
56  double Lambda = E*Nu/(1.0+Nu)/(1.0-2.0*Nu);
57  double Mu = E/2.0/(1.0+Nu);
58 
59  /// Square of the frequency of the time dependence
60  double Omega_sq = 0.5;
61 
62  /// Number of elements in r-direction
63  unsigned Nr = 5;
64 
65  /// Number of elements in z-direction
66  unsigned Nz = 10;
67 
68  /// Length of domain in r direction
69  double Lr = 1.0;
70 
71  /// Length of domain in z-direction
72  double Lz = 2.0;
73 
74  /// Set up min r coordinate
75  double Rmin = 0.1;
76 
77  /// Set up min z coordinate
78  double Zmin = 0.3;
79 
80  /// Set up max r coordinate
81  double Rmax = Rmin+Lr;
82 
83  /// Set up max z coordinate
84  double Zmax = Zmin+Lz;
85 
86  /// The traction function at r=Rmin: (t_r, t_z, t_theta)
87  void boundary_traction(const double &time,
88  const Vector<double> &x,
89  const Vector<double> &n,
90  Vector<double> &result)
91  {
92  result[0] = cos(time)*(-6.0*pow(x[0],2)*Mu*cos(x[1])-
93  Lambda*(4.0*pow(x[0],2)+pow(x[0],3))*cos(x[1]));
94  result[1] = cos(time)*(-Mu*(3.0*pow(x[0],2)-pow(x[0],3))*sin(x[1]));
95  result[2] = cos(time)*(-Mu*pow(x[0],2)*(2*pow(x[1],3)));
96  }
97 
98  /// \short The body force function; returns vector of doubles
99  /// in the order (b_r, b_z, b_theta)
100  void body_force(const double &time,
101  const Vector<double> &x,
102  Vector<double> &result)
103  {
104  result[0] = cos(time)*(
105  x[0]*(-cos(x[1])*
106  (Lambda*(8.0+3.0*x[0])-
107  Mu*(-16.0+x[0]*(x[0]-3.0))+pow(x[0],2)*Omega_sq)));
108  result[1] = cos(time)*(
109  x[0]*sin(x[1])*(Mu*(-9.0)+
110  4.0*x[0]*(Lambda+Mu)+pow(x[0],2)*
111  (Lambda+2.0*Mu-Omega_sq)));
112  result[2] = cos(time)*(
113  -x[0]*(8.0*Mu*pow(x[1],3)+pow(x[0],2)*(pow(x[1],3)*Omega_sq+6.0*Mu*x[1])));
114  } // end of body force
115 
116  /// \short Helper function - spatial components of the exact solution in a
117  /// vector. This is necessary because we need to multiply this by different
118  /// things to obtain the velocity and acceleration
119  /// 0: u_r, 1: u_z, 2: u_theta
120  void exact_solution_th(const Vector<double> &x,
121  Vector<double> &u)
122  {
123  u[0] = pow(x[0],3)*cos(x[1]);
124  u[1] = pow(x[0],3)*sin(x[1]);
125  u[2] = pow(x[0],3)*pow(x[1],3);
126  }
127 
128 
129  //Displacement, velocity and acceleration functions which are used to provide
130  //initial values to the timestepper. We must provide a separate function for
131  //each component in each case.
132 
133  /// Calculate the time dependent form of the r-component of displacement
134  double u_r(const double &time, const Vector<double> &x)
135  {
136  Vector<double> displ(3);
137  exact_solution_th(x,displ);
138  return cos(time)*displ[0];
139  } // end_of_u_r
140 
141  /// Calculate the time dependent form of the z-component of displacement
142  double u_z(const double &time, const Vector<double> &x)
143  {
144  Vector<double> displ(3);
145  exact_solution_th(x,displ);
146  return cos(time)*displ[1];
147  }
148 
149  /// Calculate the time dependent form of the theta-component of displacement
150  double u_theta(const double &time, const Vector<double> &x)
151  {
152  Vector<double> displ(3);
153  exact_solution_th(x,displ);
154  return cos(time)*displ[2];
155  }
156 
157  /// Calculate the time dependent form of the r-component of velocity
158  double d_u_r_dt(const double &time, const Vector<double> &x)
159  {
160  Vector<double> displ(3);
161  exact_solution_th(x,displ);
162  return -sin(time)*displ[0];
163  }
164 
165  /// Calculate the time dependent form of the z-component of velocity
166  double d_u_z_dt(const double &time, const Vector<double> &x)
167  {
168  Vector<double> displ(3);
169  exact_solution_th(x,displ);
170  return -sin(time)*displ[1];
171  }
172 
173  /// Calculate the time dependent form of the theta-component of velocity
174  double d_u_theta_dt(const double &time, const Vector<double> &x)
175  {
176  Vector<double> displ(3);
177  exact_solution_th(x,displ);
178  return -sin(time)*displ[2];
179  }
180 
181  /// Calculate the time dependent form of the r-component of acceleration
182  double d2_u_r_dt2(const double &time, const Vector<double> &x)
183  {
184  Vector<double> displ(3);
185  exact_solution_th(x,displ);
186  return -cos(time)*displ[0];
187  }
188 
189  /// Calculate the time dependent form of the z-component of acceleration
190  double d2_u_z_dt2(const double &time, const Vector<double> &x)
191  {
192  Vector<double> displ(3);
193  exact_solution_th(x,displ);
194  return -cos(time)*displ[1];
195  }
196 
197  /// Calculate the time dependent form of the theta-component of acceleration
198  double d2_u_theta_dt2(const double &time, const Vector<double> &x)
199  {
200  Vector<double> displ(3);
201  exact_solution_th(x,displ);
202  return -cos(time)*displ[2];
203  }
204 
205  /// The exact solution in a vector:
206  /// 0: u_r, 1: u_z, 2: u_theta and their 1st and 2nd derivs
207  void exact_solution(const double &time,
208  const Vector<double> &x,
209  Vector<double> &u)
210  {
211  u[0]=u_r(time,x);
212  u[1]=u_z(time,x);
213  u[2]=u_theta(time,x);
214 
215  u[3]=d_u_r_dt(time,x);
216  u[4]=d_u_z_dt(time,x);
217  u[5]=d_u_theta_dt(time,x);
218 
219  u[6]=d2_u_r_dt2(time,x);
220  u[7]=d2_u_z_dt2(time,x);
221  u[8]=d2_u_theta_dt2(time,x);
222  } // end_of_exact_solution
223 
224 } // end_of_Global_Parameters_namespace
225 
226 //===start_of_problem_class=============================================
227 /// Class to validate time harmonic linear elasticity (Fourier
228 /// decomposed)
229 //======================================================================
230 template<class ELEMENT, class TIMESTEPPER>
232 {
233 public:
234 
235  /// \short Constructor: Pass number of elements in r and z directions,
236  /// boundary locations and whether we are doing an impulsive start or not
238 
239  /// Update before solve is empty
241 
242  /// Update after solve is empty
244 
245  /// Actions before implicit timestep
247  {
248  // Just need to update the boundary conditions
249  set_boundary_conditions();
250  }
251 
252  /// \short Set the initial conditions, either for an impulsive start or
253  /// with history values for the time stepper
254  void set_initial_conditions();
255 
256  /// Set the boundary conditions
257  void set_boundary_conditions();
258 
259  /// Doc the solution
260  void doc_solution(DocInfo& doc_info);
261 
262 private:
263 
264  /// Allocate traction elements on the bottom surface
265  void assign_traction_elements();
266 
267  /// Pointer to the bulk mesh
269 
270  /// Pointer to the mesh of traction elements
272 }; // end_of_problem_class
273 
274 
275 //===start_of_constructor=============================================
276 /// Problem constructor: Pass number of elements in coordinate
277 /// directions and size of domain.
278 //====================================================================
279 template<class ELEMENT, class TIMESTEPPER>
282 {
283  //Allocate the timestepper
284  add_time_stepper_pt(new TIMESTEPPER());
285 
286  //Now create the mesh
287  Bulk_mesh_pt = new RectangularQuadMesh<ELEMENT>(
294  time_stepper_pt());
295 
296  //Create the surface mesh of traction elements
297  Surface_mesh_pt=new Mesh;
298  assign_traction_elements();
299 
300  //Set the boundary conditions
301  set_boundary_conditions();
302 
303  // Complete the problem setup to make the elements fully functional
304 
305  // Loop over the elements
306  unsigned n_el = Bulk_mesh_pt->nelement();
307  for(unsigned e=0;e<n_el;e++)
308  {
309  // Cast to a bulk element
310  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
311 
312  // Set the body force
313  el_pt->body_force_fct_pt() = &Global_Parameters::body_force;
314 
315  // Set the pointer to Poisson's ratio
316  el_pt->nu_pt() = &Global_Parameters::Nu;
317 
318  // Set the pointer to non-dim Young's modulus
319  el_pt->youngs_modulus_pt() = &Global_Parameters::E;
320 
321  // Set the pointer to the Lambda parameter
322  el_pt->lambda_sq_pt() = &Global_Parameters::Omega_sq;
323 
324  }// end_loop_over_elements
325 
326  // Loop over the traction elements
327  unsigned n_traction = Surface_mesh_pt->nelement();
328  for(unsigned e=0;e<n_traction;e++)
329  {
330  // Cast to a surface element
331  AxisymmetricLinearElasticityTractionElement<ELEMENT>*
332  el_pt =
333  dynamic_cast<AxisymmetricLinearElasticityTractionElement
334  <ELEMENT>* >(Surface_mesh_pt->element_pt(e));
335 
336  // Set the applied traction
337  el_pt->traction_fct_pt() = &Global_Parameters::boundary_traction;
338 
339  }// end_loop_over_traction_elements
340 
341  // Add the submeshes to the problem
342  add_sub_mesh(Bulk_mesh_pt);
343  add_sub_mesh(Surface_mesh_pt);
344 
345  // Now build the global mesh
346  build_global_mesh();
347 
348  // Assign equation numbers
349  cout << assign_eqn_numbers() << " equations assigned" << std::endl;
350 
351 
352  //{
353  // // Create animation
354  // unsigned ntime=20;
355  // double displ_r=0;
356  // char filename[30];
357  // unsigned n_node=Bulk_mesh_pt->nnode();
358  // Node* nod_pt=0;
359 
360  // for (unsigned it=0;it<ntime;it++)
361  // {
362  // sprintf(filename,"animation%i.dat",it);
363  // Global_Parameters::Output_stream.open(filename);
364  // double t=(2*MathematicalConstants::Pi)*(double(it)/(ntime-1));
365  // std::cout << t << std::endl;
366  //
367  // // Loop over nodes of the mesh
368  // for(unsigned j=0;j<n_node;j++)
369  // {
370  // nod_pt=Bulk_mesh_pt->node_pt(j);
371  // Vector<double> x(2);
372  // x[0]=nod_pt->x(0);
373  // x[1]=nod_pt->x(1);
374 
375  // displ_r = Global_Parameters::u_r(t,x);
376  // Global_Parameters::Output_stream << x[0] << ' ';
377  // Global_Parameters::Output_stream << x[1] << ' ';
378  // Global_Parameters::Output_stream << displ_r << std::endl;
379  // }
380  // Global_Parameters::Output_stream.close();
381  // }
382  //} //end_of_animation
383 
384 } // end_of_constructor
385 
386 
387 //===start_of_traction===============================================
388 /// Make traction elements along the boundary r=Rmin
389 //===================================================================
390 template<class ELEMENT, class TIMESTEPPER>
393 {
394  unsigned bound, n_neigh;
395 
396  // How many bulk elements are next to boundary 3
397  bound=3;
398  n_neigh = Bulk_mesh_pt->nboundary_element(bound);
399 
400  // Now loop over bulk elements and create the face elements
401  for(unsigned n=0;n<n_neigh;n++)
402  {
403  // Create the face element
404  FiniteElement *traction_element_pt
405  = new AxisymmetricLinearElasticityTractionElement<ELEMENT>
406  (Bulk_mesh_pt->boundary_element_pt(bound,n),
407  Bulk_mesh_pt->face_index_at_boundary(bound,n));
408 
409  // Add to mesh
410  Surface_mesh_pt->add_element_pt(traction_element_pt);
411  }
412 
413 } // end of assign_traction_elements
414 
415 //===start_of_set_initial_conditions=================================
416 /// Set the initial conditions (history values)
417 //===================================================================
418 template<class ELEMENT, class TIMESTEPPER>
421 {
422  // Upcast the timestepper to the specific type we have
423  TIMESTEPPER* timestepper_pt =
424  dynamic_cast<TIMESTEPPER*>(time_stepper_pt());
425 
426  // By default do a non-impulsive start and provide initial conditions
427  bool impulsive_start=false;
428 
429  if(impulsive_start)
430  {
431  // Number of nodes in the bulk mesh
432  unsigned n_node = Bulk_mesh_pt->nnode();
433 
434  // Loop over all nodes in the bulk mesh
435  for(unsigned inod=0;inod<n_node;inod++)
436  {
437  // Pointer to node
438  Node* nod_pt = Bulk_mesh_pt->node_pt(inod);
439 
440  // Get nodal coordinates
441  Vector<double> x(2);
442  x[0] = nod_pt->x(0);
443  x[1] = nod_pt->x(1);
444 
445  // Assign zero solution at t=0
446  nod_pt->set_value(0,0);
447  nod_pt->set_value(1,0);
448  nod_pt->set_value(2,0);
449 
450  // Set the impulsive initial values in the timestepper
451  timestepper_pt->assign_initial_values_impulsive(nod_pt);
452  }
453  } // end_of_impulsive_start
454  else // Smooth start
455  {
456  // Storage for pointers to the functions defining the displacement,
457  // velocity and acceleration components
458  Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
459  initial_value_fct(3);
460  Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
461  initial_veloc_fct(3);
462  Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
463  initial_accel_fct(3);
464 
465  // Set the displacement function pointers
466  initial_value_fct[0]=&Global_Parameters::u_r;
467  initial_value_fct[1]=&Global_Parameters::u_z;
468  initial_value_fct[2]=&Global_Parameters::u_theta;
469 
470  // Set the velocity function pointers
471  initial_veloc_fct[0]=&Global_Parameters::d_u_r_dt;
472  initial_veloc_fct[1]=&Global_Parameters::d_u_z_dt;
473  initial_veloc_fct[2]=&Global_Parameters::d_u_theta_dt;
474 
475  // Set the acceleration function pointers
476  initial_accel_fct[0]=&Global_Parameters::d2_u_r_dt2;
477  initial_accel_fct[1]=&Global_Parameters::d2_u_z_dt2;
478  initial_accel_fct[2]=&Global_Parameters::d2_u_theta_dt2;
479 
480  // Number of nodes in the bulk mesh
481  unsigned n_node = Bulk_mesh_pt->nnode();
482 
483  // Loop over all nodes in bulk mesh
484  for(unsigned inod=0;inod<n_node;inod++)
485  {
486  // Pointer to node
487  Node* nod_pt = Bulk_mesh_pt->node_pt(inod);
488 
489  // Assign the history values
490  timestepper_pt->assign_initial_data_values(nod_pt,
491  initial_value_fct,
492  initial_veloc_fct,
493  initial_accel_fct);
494  } // end_of_loop_over_nodes
495 
496  // Paranoid checking of history values
497  double err_max=0.0;
498 
499  // Loop over all nodes in bulk mesh
500  for(unsigned jnod=0;jnod<n_node;jnod++)
501  {
502  // Pointer to node
503  Node* nod_pt=Bulk_mesh_pt->node_pt(jnod);
504 
505  // Get nodal coordinates
506  Vector<double> x(2);
507  x[0]=nod_pt->x(0);
508  x[1]=nod_pt->x(1);
509 
510  // Get exact displacements
511  double u_r_exact=
512  Global_Parameters::u_r(time_pt()->time(),x);
513  double u_z_exact=
514  Global_Parameters::u_z(time_pt()->time(),x);
515  double u_theta_exact=
516  Global_Parameters::u_theta(time_pt()->time(),x);
517 
518  // Get exact velocities
519  double d_u_r_dt_exact=
520  Global_Parameters::d_u_r_dt(time_pt()->time(),x);
521  double d_u_z_dt_exact=
522  Global_Parameters::d_u_z_dt(time_pt()->time(),x);
523  double d_u_theta_dt_exact=
524  Global_Parameters::d_u_theta_dt(time_pt()->time(),x);
525 
526  // Get exact accelerations
527  double d2_u_r_dt2_exact=
528  Global_Parameters::d2_u_r_dt2(time_pt()->time(),x);
529  double d2_u_z_dt2_exact=
530  Global_Parameters::d2_u_z_dt2(time_pt()->time(),x);
531  double d2_u_theta_dt2_exact=
532  Global_Parameters::d2_u_theta_dt2(time_pt()->time(),x);
533 
534  // Get Newmark approximations for:
535  // zero-th time derivatives
536  double u_r_fe=timestepper_pt->time_derivative(0,nod_pt,0);
537  double u_z_fe=timestepper_pt->time_derivative(0,nod_pt,1);
538  double u_theta_fe=timestepper_pt->time_derivative(0,nod_pt,2);
539 
540  // first time derivatives
541  double d_u_r_dt_fe=timestepper_pt->time_derivative(1,nod_pt,0);
542  double d_u_z_dt_fe=timestepper_pt->time_derivative(1,nod_pt,1);
543  double d_u_theta_dt_fe=timestepper_pt->time_derivative(1,nod_pt,2);
544 
545  // second time derivatives
546  double d2_u_r_dt2_fe=timestepper_pt->time_derivative(2,nod_pt,0);
547  double d2_u_z_dt2_fe=timestepper_pt->time_derivative(2,nod_pt,1);
548  double d2_u_theta_dt2_fe=timestepper_pt->time_derivative(2,nod_pt,2);
549 
550  // Calculate the error as the norm of all the differences between the
551  // Newmark approximations and the 'exact' (numerical) expressions
552  double error=sqrt(pow(u_r_exact-u_r_fe,2)+
553  pow(u_z_exact-u_z_fe,2)+
554  pow(u_theta_exact-u_theta_fe,2)+
555  pow(d_u_r_dt_exact-d_u_r_dt_fe,2)+
556  pow(d_u_z_dt_exact-d_u_z_dt_fe,2)+
557  pow(d_u_theta_dt_exact-d_u_theta_dt_fe,2)+
558  pow(d2_u_r_dt2_exact-d2_u_r_dt2_fe,2)+
559  pow(d2_u_z_dt2_exact-d2_u_z_dt2_fe,2)+
560  pow(d2_u_theta_dt2_exact-d2_u_theta_dt2_fe,2));
561 
562  // If there is an error greater than one previously seen, keep hold of it
563  if(error>err_max)
564  {
565  err_max=error;
566  }
567  } // end of loop over nodes
568  std::cout << "Max error in assignment of initial conditions "
569  << err_max << std::endl;
570  }
571 } // end_of_set_initial_conditions
572 
573 //==start_of_set_boundary_conditions======================================
574 /// Set the boundary conditions
575 //========================================================================
576 template<class ELEMENT, class TIMESTEPPER>
579 {
580  // Set the boundary conditions for this problem: All nodes are
581  // free by default -- just pin & set the ones that have Dirichlet
582  // conditions here
583 
584  // storage for nodal position
585  Vector<double> x(2);
586 
587  // Storage for prescribed displacements
588  Vector<double> u(9);
589 
590  // Storage for pointers to the functions defining the displacement,
591  // velocity and acceleration components
592  Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
593  initial_value_fct(3);
594  Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
595  initial_veloc_fct(3);
596  Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
597  initial_accel_fct(3);
598 
599  // Set the displacement function pointers
600  initial_value_fct[0]=&Global_Parameters::u_r;
601  initial_value_fct[1]=&Global_Parameters::u_z;
602  initial_value_fct[2]=&Global_Parameters::u_theta;
603 
604  // Set the velocity function pointers
605  initial_veloc_fct[0]=&Global_Parameters::d_u_r_dt;
606  initial_veloc_fct[1]=&Global_Parameters::d_u_z_dt;
607  initial_veloc_fct[2]=&Global_Parameters::d_u_theta_dt;
608 
609  // Set the acceleration function pointers
610  initial_accel_fct[0]=&Global_Parameters::d2_u_r_dt2;
611  initial_accel_fct[1]=&Global_Parameters::d2_u_z_dt2;
612  initial_accel_fct[2]=&Global_Parameters::d2_u_theta_dt2;
613 
614 
615  // Now set displacements on boundaries 0 (z=Zmin),
616  //------------------------------------------------
617  // 1 (r=Rmax) and 2 (z=Zmax)
618  //--------------------------
619  for (unsigned ibound=0;ibound<=2;ibound++)
620  {
621  unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
622  for (unsigned inod=0;inod<num_nod;inod++)
623  {
624  // Get pointer to node
625  Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(ibound,inod);
626 
627  // Pinned in r, z and theta
628  nod_pt->pin(0);nod_pt->pin(1);nod_pt->pin(2);
629 
630  // Direct assigment of just the current values...
631  bool use_direct_assigment=true;
632  if (use_direct_assigment)
633  {
634  // get r and z coordinates
635  x[0]=nod_pt->x(0);
636  x[1]=nod_pt->x(1);
637 
638  // Compute the value of the exact solution at the nodal point
639  Global_Parameters::exact_solution(time_pt()->time(),x,u);
640 
641  // Set the displacements
642  nod_pt->set_value(0,u[0]);
643  nod_pt->set_value(1,u[1]);
644  nod_pt->set_value(2,u[2]);
645  }
646  // ...or the history values too:
647  else
648  {
649  // Upcast the timestepper to the specific type we have
650  TIMESTEPPER* timestepper_pt =
651  dynamic_cast<TIMESTEPPER*>(time_stepper_pt());
652 
653  // Assign the history values
654  timestepper_pt->assign_initial_data_values(nod_pt,
655  initial_value_fct,
656  initial_veloc_fct,
657  initial_accel_fct);
658  }
659  }
660  } // end_of_loop_over_boundary_nodes
661 } // end_of_set_boundary_conditions
662 
663 //==start_of_doc_solution=================================================
664 /// Doc the solution
665 //========================================================================
666 template<class ELEMENT, class TIMESTEPPER>
668 doc_solution(DocInfo& doc_info)
669 {
670  ofstream some_file;
671  char filename[100];
672 
673  // Number of plot points
674  unsigned npts=10;
675 
676  // Output solution
677  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
678  doc_info.number());
679  some_file.open(filename);
680  Bulk_mesh_pt->output(some_file,npts);
681  some_file.close();
682 
683  // Output exact solution
684  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
685  doc_info.number());
686  some_file.open(filename);
687  Bulk_mesh_pt->output_fct(some_file,npts,time_pt()->time(),
689  some_file.close();
690 
691  // Doc error
692  double error=0.0;
693  double norm=0.0;
694  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
695  doc_info.number());
696  some_file.open(filename);
697  Bulk_mesh_pt->compute_error(some_file,
699  time_pt()->time(),
700  error,norm);
701  some_file.close();
702 
703  // Doc error norm:
704  cout << "\nNorm of error: " << sqrt(error) << std::endl;
705  cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
706  cout << std::endl;
707 } // end_of_doc_solution
708 
709 //======================================================================
710 // End of Axisymmetric problem class
711 //======================================================================
712 
713 
714 
715 //===start_of_main======================================================
716 /// Driver code
717 //======================================================================
718 int main(int argc, char* argv[])
719 {
720  // Store command line arguments
721  CommandLineArgs::setup(argc,argv);
722 
723  // Define possible command line arguments and parse the ones that
724  // were actually specified
725 
726  // Validation?
727  CommandLineArgs::specify_command_line_flag("--validation");
728 
729  // Parse command line
730  CommandLineArgs::parse_and_assign();
731 
732  // Doc what has actually been specified on the command line
733  CommandLineArgs::doc_specified_flags();
734 
735  // Set up doc info
736  DocInfo doc_info;
737 
738  // Set output directory
739  doc_info.set_directory("RESLT");
740 
741  // Time dependent problem instance
743  <QAxisymmetricLinearElasticityElement<3>, Newmark<1> > problem;
744 
745  // Set the initial time to t=0
746  problem.time_pt()->time()=0.0;
747 
748  // Set and initialise timestep
749  double dt=0;
750 
751  // If we're validating, use a larger timestep so that we can do fewer steps
752  // before reaching interesting behaviour
753  if(CommandLineArgs::command_line_flag_has_been_set("--validation"))
754  {
755  dt=0.1;
756  }
757  else // Otherwise use a small timestep
758  {
759  dt=0.01;
760  }
761  problem.time_pt()->initialise_dt(dt);
762 
763  // Set the initial conditions
764  problem.set_initial_conditions();
765 
766  // Doc the initial conditions and increment the doc_info number
767  problem.doc_solution(doc_info);
768  doc_info.number()++;
769 
770  // Find the number of timesteps to perform
771  unsigned nstep=0;
772 
773  // If we're validating, only do a few timesteps; otherwise do a whole period
774  if(CommandLineArgs::command_line_flag_has_been_set("--validation"))
775  {
776  nstep=5;
777  }
778  else // Otherwise calculate based on timestep
779  {
780  // Solve for one full period
781  double t_max=2*MathematicalConstants::Pi;
782 
783  nstep=unsigned(t_max/dt);
784  } //end_of_calculate_number_of_timesteps
785 
786  // Do the timestepping
787  for(unsigned istep=0;istep<nstep;istep++)
788  {
789  // Solve for this timestep
790  problem.unsteady_newton_solve(dt);
791 
792  // Doc the solution and increment doc_info number
793  problem.doc_solution(doc_info);
794  doc_info.number()++;
795  }
796 } // end_of_main
797 
double Omega_sq
Square of the frequency of the time dependence.
Definition: cylinder.cc:60
double d_u_theta_dt(const double &time, const Vector< double > &x)
Calculate the time dependent form of the theta-component of velocity.
Definition: cylinder.cc:174
double d2_u_r_dt2(const double &time, const Vector< double > &x)
Calculate the time dependent form of the r-component of acceleration.
Definition: cylinder.cc:182
double d_u_r_dt(const double &time, const Vector< double > &x)
Calculate the time dependent form of the r-component of velocity.
Definition: cylinder.cc:158
void body_force(const double &time, const Vector< double > &x, Vector< double > &result)
The body force function; returns vector of doubles in the order (b_r, b_z, b_theta) ...
Definition: cylinder.cc:100
double u_z(const double &time, const Vector< double > &x)
Calculate the time dependent form of the z-component of displacement.
Definition: cylinder.cc:142
void boundary_traction(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
The traction function at r=Rmin: (t_r, t_z, t_theta)
Definition: cylinder.cc:87
AxisymmetricLinearElasticityProblem()
Constructor: Pass number of elements in r and z directions, boundary locations and whether we are doi...
Definition: cylinder.cc:281
double Rmin
Set up min r coordinate.
Definition: cylinder.cc:75
unsigned Nr
Number of elements in r-direction.
Definition: cylinder.cc:63
void exact_solution_th(const Vector< double > &x, Vector< double > &u)
Helper function - spatial components of the exact solution in a vector. This is necessary because we ...
Definition: cylinder.cc:120
double d2_u_z_dt2(const double &time, const Vector< double > &x)
Calculate the time dependent form of the z-component of acceleration.
Definition: cylinder.cc:190
double E
Define the non-dimensional Young's modulus.
Definition: cylinder.cc:53
void exact_solution(const double &time, const Vector< double > &x, Vector< double > &u)
Definition: cylinder.cc:207
void assign_traction_elements()
Allocate traction elements on the bottom surface.
Definition: cylinder.cc:392
double d2_u_theta_dt2(const double &time, const Vector< double > &x)
Calculate the time dependent form of the theta-component of acceleration.
Definition: cylinder.cc:198
void set_initial_conditions()
Set the initial conditions, either for an impulsive start or with history values for the time stepper...
Definition: cylinder.cc:420
void actions_after_newton_solve()
Update after solve is empty.
Definition: cylinder.cc:243
void actions_before_newton_solve()
Update before solve is empty.
Definition: cylinder.cc:240
unsigned Nz
Number of elements in z-direction.
Definition: cylinder.cc:66
double Lr
Length of domain in r direction.
Definition: cylinder.cc:69
double Zmax
Set up max z coordinate.
Definition: cylinder.cc:84
void set_boundary_conditions()
Set the boundary conditions.
Definition: cylinder.cc:578
Mesh * Surface_mesh_pt
Pointer to the mesh of traction elements.
Definition: cylinder.cc:271
double u_r(const double &time, const Vector< double > &x)
Calculate the time dependent form of the r-component of displacement.
Definition: cylinder.cc:134
double Lz
Length of domain in z-direction.
Definition: cylinder.cc:72
double Nu
Define Poisson's ratio Nu.
Definition: cylinder.cc:50
double Lambda
Lame parameters.
Definition: cylinder.cc:56
void actions_before_implicit_timestep()
Actions before implicit timestep.
Definition: cylinder.cc:246
double Rmax
Set up max r coordinate.
Definition: cylinder.cc:81
double d_u_z_dt(const double &time, const Vector< double > &x)
Calculate the time dependent form of the z-component of velocity.
Definition: cylinder.cc:166
void doc_solution(DocInfo &doc_info)
Doc the solution.
Definition: cylinder.cc:668
int main(int argc, char *argv[])
Driver code.
Definition: cylinder.cc:718
Mesh * Bulk_mesh_pt
Pointer to the bulk mesh.
Definition: cylinder.cc:268
double u_theta(const double &time, const Vector< double > &x)
Calculate the time dependent form of the theta-component of displacement.
Definition: cylinder.cc:150
double Zmin
Set up min z coordinate.
Definition: cylinder.cc:78