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