disk_oscillation.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 elasticity problem
31 
32 // oomph-lib includes
33 #include "generic.h"
34 #include "solid.h"
35 #include "oomph_crbond_bessel.h"
36 
37 //Need to instantiate templated mesh
38 #include "meshes/quarter_circle_sector_mesh.h"
39 
40 using namespace std;
41 
42 using namespace oomph;
43 
44 
45 ///////////////////////////////////////////////////////////////////////
46 ///////////////////////////////////////////////////////////////////////
47 ///////////////////////////////////////////////////////////////////////
48 
49 
50 
51 //=======start_namespace==========================================
52 /// Global variables
53 //================================================================
54 namespace Global_Physical_Variables
55 {
56  /// Poisson's ratio
57  double Nu=0.3;
58 
59  /// Timescale ratio
60  double Lambda_sq=(1.0-Nu)/((1.0+Nu)*(1.0-2.0*Nu));
61 
62  /// Pointer to constitutive law
63  ConstitutiveLaw* Constitutive_law_pt=0;
64 
65  /// \short Multiplier for inertia terms (needed for consistent assignment
66  /// of initial conditions in Newmark scheme)
67  double multiplier(const Vector<double>& xi)
68  {
70  }
71 
72 
73 } // end namespace
74 
75 
76 ///////////////////////////////////////////////////////////////////////
77 ///////////////////////////////////////////////////////////////////////
78 // Axisymmetrically oscillating disk
79 ///////////////////////////////////////////////////////////////////////
80 ///////////////////////////////////////////////////////////////////////
81 
82 
83 
84 //================disk_as_geom_object======================================
85 /// \short Axisymmetrially oscillating disk with displacement
86 /// field according to linear elasticity.
87 //=========================================================================
88 class AxisymOscillatingDisk : public GeomObject
89 {
90 
91 public:
92 
93  /// \short Constructor: 2 Lagrangian coordinate, 2 Eulerian coords. Pass
94  /// amplitude of oscillation, Poisson ratio nu, and pointer to
95  /// global timestepper.
96  AxisymOscillatingDisk(const double& ampl, const double& nu,
97  TimeStepper* time_stepper_pt);
98 
99  /// Destructor (empty)
101 
102  /// \short Position vector at Lagrangian coordinate xi at present
103  /// time
104  void position(const Vector<double>& xi, Vector<double>& r) const;
105 
106  ///\short Parametrised velocity on object at current time: veloc = d r(xi)/dt.
107  void veloc(const Vector<double>& xi, Vector<double>& veloc);
108 
109  /// \short Parametrised acceleration on object at current time:
110  /// accel = d^2 r(xi)/dt^2.
111  void accel(const Vector<double>& xi, Vector<double>& accel);
112 
113  /// \short Parametrised j-th time-derivative on object at current time:
114  /// \f$ \frac{d^{j} r(\zeta)}{dt^j} \f$.
115  void dposition_dt(const Vector<double>& xi, const unsigned& j,
116  Vector<double>& drdt)
117  {
118  switch (j)
119  {
120  // Current position
121  case 0:
122  position(xi,drdt);
123  break;
124 
125  // Velocity:
126  case 1:
127  veloc(xi,drdt);
128  break;
129 
130  // Acceleration:
131  case 2:
132  accel(xi,drdt);
133  break;
134 
135  default:
136  std::ostringstream error_message;
137  error_message << j << "th derivative not implemented\n";
138 
139  throw OomphLibError(error_message.str(),
140  OOMPH_CURRENT_FUNCTION,
141  OOMPH_EXCEPTION_LOCATION);
142  }
143  }
144 
145 
146  /// \short Residual of dispersion relation for use in black-box Newton method
147  /// which requires global (or static) functions.
148  /// Poisson's ratio is passed as parameter.
149  static void residual_for_dispersion(const Vector<double>& param,
150  const Vector<double>& omega,
151  Vector<double>& residual);
152 
153 private:
154 
155  /// Amplitude of oscillation
156  double Ampl;
157 
158  /// Period of oscillation
159  double T;
160 
161  /// Poisson ratio nu
162  double Nu;
163 
164  /// Eigenfrequency
165  double Omega;
166 
167 }; // end disk_as_geom_object
168 
169 
170 
171 
172 //==============ic_constructor============================================
173 /// Constructor: 2 Lagrangian coordinates, 2 Eulerian coords. Pass
174 /// amplitude of oscillation, Poisson ratio nu, and pointer to
175 /// global timestepper.
176 //========================================================================
178  const double& nu,
179  TimeStepper* time_stepper_pt) :
180  GeomObject(2,2,time_stepper_pt), Ampl(ampl), Nu(nu)
181 {
182  // Parameters for dispersion relation
183  Vector<double> param(1);
184  param[0]=Nu;
185 
186  // Initial guess for eigenfrequency
187  Vector<double> omega(1);
188  omega[0]=2.0;
189 
190  // Find eigenfrequency from black box Newton solver
191  BlackBoxFDNewtonSolver::black_box_fd_newton_solve(residual_for_dispersion,
192  param,omega);
193 
194  // Assign eigenfrequency
195  Omega=omega[0];
196 
197  // Assign/doc period of oscillation
198  T=2.0*MathematicalConstants::Pi/Omega;
199 
200  std::cout << "Period of oscillation: " << T << std::endl;
201 
202 }
203 
204 
205 
206 //===============start_position===========================================
207 /// \short Position Vector at Lagrangian coordinate xi at present
208 /// time
209 //========================================================================
210 void AxisymOscillatingDisk::position(const Vector<double>& xi,
211  Vector<double>& r) const
212 {
213  // Parameter values at present time
214  double time=Geom_object_time_stepper_pt->time_pt()->time();
215 
216  // Radius in Lagrangian coordinates
217  double lagr_radius=sqrt( pow(xi[0],2) + pow(xi[1],2) );
218 
219  if (lagr_radius<1.0e-12)
220  {
221  // Position Vector
222  r[0]=0.0;
223  r[1]=0.0;
224  }
225  else
226  {
227  // Bessel fcts J_0(x), J_1(x), Y_0(x), Y_1(x) and their derivatives
228  double j0,j1,y0,y1,j0p,j1p,y0p,y1p;
229  CRBond_Bessel::bessjy01a(Omega*lagr_radius,j0,j1,y0,y1,j0p,j1p,y0p,y1p);
230 
231  // Displacement field
232  double u=Ampl*j1*sin(2.0*MathematicalConstants::Pi*time/T);
233 
234  // Position Vector
235  r[0]=(xi[0]+xi[0]/lagr_radius*u);
236  r[1]=(xi[1]+xi[1]/lagr_radius*u);
237  }
238 
239 } //end position
240 
241 
242 
243 
244 //========================================================================
245 ///\short Parametrised velocity on object at current time:
246 /// veloc = d r(xi)/dt.
247 //========================================================================
248 void AxisymOscillatingDisk::veloc(const Vector<double>& xi,
249  Vector<double>& veloc)
250 {
251  // Parameter values at present time
252  double time=Geom_object_time_stepper_pt->time_pt()->time();
253 
254  // Radius in Lagrangian coordinates
255  double lagr_radius=sqrt(pow(xi[0],2)+pow(xi[1],2));
256 
257  if (lagr_radius<1.0e-12)
258  {
259  // Veloc vector
260  veloc[0]=0.0;
261  veloc[1]=0.0;
262  }
263  else
264  {
265  // Bessel fcts J_0(x), J_1(x), Y_0(x), Y_1(x) and their derivatives
266  double j0,j1,y0,y1,j0p,j1p,y0p,y1p;
267  CRBond_Bessel::bessjy01a(Omega*lagr_radius,j0,j1,y0,y1,j0p,j1p,y0p,y1p);
268 
269  // Deriv of displacement field
270  double udot=2.0*MathematicalConstants::Pi/T*Ampl*j1*
271  cos(2.0*MathematicalConstants::Pi*time/T);
272 
273  // Veloc
274  veloc[0]=(xi[0]/lagr_radius*udot);
275  veloc[1]=(xi[1]/lagr_radius*udot);
276  }
277 }
278 
279 
280 
281 
282 
283 //========================================================================
284 /// Parametrised acceleration on object at current time:
285 /// accel = d^2 r(xi)/dt^2.
286 //========================================================================
287 void AxisymOscillatingDisk::accel(const Vector<double>& xi,
288  Vector<double>& accel)
289 {
290  // Parameter values at present time
291  double time=Geom_object_time_stepper_pt->time_pt()->time();
292 
293  // Radius in Lagrangian coordinates
294  double lagr_radius=sqrt(pow(xi[0],2)+pow(xi[1],2));
295 
296 
297  if (lagr_radius<1.0e-12)
298  {
299  // Veloc vector
300  accel[0]=0.0;
301  accel[1]=0.0;
302  }
303  else
304  {
305  // Bessel fcts J_0(x), J_1(x), Y_0(x), Y_1(x) and their derivatives
306  double j0,j1,y0,y1,j0p,j1p,y0p,y1p;
307  CRBond_Bessel::bessjy01a(Omega*lagr_radius,j0,j1,y0,y1,j0p,j1p,y0p,y1p);
308 
309  // Deriv of displacement field
310  double udotdot=-pow(2.0*MathematicalConstants::Pi/T,2)*Ampl*j1*
311  sin(2.0*MathematicalConstants::Pi*time/T);
312 
313  // Veloc
314  accel[0]=(xi[0]/lagr_radius*udotdot);
315  accel[1]=(xi[1]/lagr_radius*udotdot);
316 
317  }
318 }
319 
320 
321 
322 //======================start_of_dispersion===============================
323 /// Residual of dispersion relation for use in black box Newton method
324 /// which requires global (or static) functions.
325 /// Poisson's ratio is passed as parameter.
326 //========================================================================
328  const Vector<double>& param, const Vector<double>& omega,
329  Vector<double>& residual)
330 {
331  // Extract parameters
332  double nu=param[0];
333 
334  // Argument of various Bessel functions
335  double arg=omega[0];
336 
337  // Bessel fcts J_0(x), J_1(x), Y_0(x), Y_1(x) and their derivatives
338  double j0,j1,y0,y1,j0p,j1p,y0p,y1p;
339  CRBond_Bessel::bessjy01a(arg,j0,j1,y0,y1,j0p,j1p,y0p,y1p);
340 
341  // Residual of dispersion relation
342  residual[0]=nu/(1.0-2.0*nu)*(j1+(j0-j1/omega[0])*omega[0])+
343  (j0-j1/omega[0])*omega[0];
344 
345 }
346 
347 
348 
349 ///////////////////////////////////////////////////////////////////////
350 ///////////////////////////////////////////////////////////////////////
351 ///////////////////////////////////////////////////////////////////////
352 
353 
354 
355 //======================start_mesh================================
356 /// Elastic quarter circle sector mesh: We "upgrade"
357 /// the RefineableQuarterCircleSectorMesh to become an
358 /// SolidMesh and equate the Eulerian and Lagrangian coordinates,
359 /// thus making the domain represented by the mesh the stress-free
360 /// configuration.
361 //================================================================
362 template <class ELEMENT>
364  public virtual RefineableQuarterCircleSectorMesh<ELEMENT>,
365  public virtual SolidMesh
366 {
367 
368 
369 public:
370 
371  /// \short Constructor: Build mesh and copy Eulerian coords to Lagrangian
372  /// ones so that the initial configuration is the stress-free one.
374  const double& xi_lo,
375  const double& fract_mid,
376  const double& xi_hi,
377  TimeStepper* time_stepper_pt=
378  &Mesh::Default_TimeStepper) :
379  RefineableQuarterCircleSectorMesh<ELEMENT>(wall_pt,xi_lo,fract_mid,xi_hi,
380  time_stepper_pt)
381  {
382 #ifdef PARANOID
383  /// Check that the element type is derived from the SolidFiniteElement
384  SolidFiniteElement* el_pt=dynamic_cast<SolidFiniteElement*>
385  (finite_element_pt(0));
386  if (el_pt==0)
387  {
388  throw OomphLibError(
389  "Element needs to be derived from SolidFiniteElement\n",
390  OOMPH_CURRENT_FUNCTION,
391  OOMPH_EXCEPTION_LOCATION);
392  }
393 #endif
394 
395  // Make the current configuration the undeformed one by
396  // setting the nodal Lagrangian coordinates to their current
397  // Eulerian ones
398  set_lagrangian_nodal_coordinates();
399  }
400 
401 };
402 
403 
404 
405 
406 
407 /////////////////////////////////////////////////////////////////////////
408 /////////////////////////////////////////////////////////////////////////
409 /////////////////////////////////////////////////////////////////////////
410 
411 
412 
413 //========start_class===================================================
414 /// \short Problem class to simulate small-amplitude oscillations of
415 /// a circular disk.
416 //======================================================================
417 template<class ELEMENT>
418 class DiskOscillationProblem : public Problem
419 {
420 
421 public:
422 
423  /// Constructor
425 
426  /// Update function (empty)
428 
429  /// Update function (empty)
431 
432  /// Access function for the solid mesh
434  {
436  Problem::mesh_pt());
437  }
438 
439  /// Run the problem: Pass number of timesteps to be performed.
440  void run(const unsigned& nstep);
441 
442  /// Doc the solution
443  void doc_solution(DocInfo& doc_info);
444 
445 private:
446 
447  /// Trace file
448  ofstream Trace_file;
449 
450  /// Vector of pointers to nodes whose position we're tracing
451  Vector<Node*> Trace_node_pt;
452 
453  /// Geometric object that specifies the initial conditions
455 
456 }; // end class
457 
458 
459 
460 
461 
462 //============start_constructor=========================================
463 /// Constructor
464 //======================================================================
465 template<class ELEMENT>
467 {
468 
469  // Allocate the timestepper: The classical Newmark scheme with
470  // two history values.
471  add_time_stepper_pt(new Newmark<2>);
472 
473 
474 
475  // GeomObject that specifies the curvilinear boundary of the
476  // circular disk
477  GeomObject* curved_boundary_pt=new Ellipse(1.0,1.0);
478 
479  //The start and end intrinsic coordinates on the geometric object
480  // that defines the curvilinear boundary of the disk
481  double xi_lo=0.0;
482  double xi_hi=2.0*atan(1.0);
483 
484  // Fraction along geometric object at which the radial dividing line
485  // is placed
486  double fract_mid=0.5;
487 
488  //Now create the mesh
490  curved_boundary_pt,xi_lo,fract_mid,xi_hi,time_stepper_pt());
491 
492 
493  // Setup trace nodes as the nodes on boundaries 0 (= horizontal symmetry
494  // boundary) and 1 (=curved boundary)
495  unsigned nnod0=mesh_pt()->nboundary_node(0);
496  unsigned nnod1=mesh_pt()->nboundary_node(1);
497  Trace_node_pt.resize(nnod0+nnod1);
498  for (unsigned j=0;j<nnod0;j++)
499  {
500 
501  Trace_node_pt[j]=mesh_pt()->boundary_node_pt(0,j);
502 
503  }
504  for (unsigned j=0;j<nnod1;j++)
505  {
506 
507  Trace_node_pt[j+nnod0]=mesh_pt()->boundary_node_pt(1,j);
508 
509  } //done choosing trace nodes
510 
511 
512  // Pin the horizontal boundary in the vertical direction
513  unsigned n_hor = mesh_pt()->nboundary_node(0);
514  for(unsigned i=0;i<n_hor;i++)
515  {
516  mesh_pt()->boundary_node_pt(0,i)->pin_position(1);
517  }
518 
519  // Pin the vertical boundary in the horizontal direction
520  unsigned n_vert = mesh_pt()->nboundary_node(2);
521  for(unsigned i=0;i<n_vert;i++)
522  {
523 
524  mesh_pt()->boundary_node_pt(2,i)->pin_position(0);
525 
526  } // done bcs
527 
528 
529  //Finish build of elements
530  unsigned n_element =mesh_pt()->nelement();
531  for(unsigned i=0;i<n_element;i++)
532  {
533  //Cast to a solid element
534  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
535 
536  //Set the constitutive law
537  el_pt->constitutive_law_pt() =
539 
540  // Set the timescale ratio
541  el_pt->lambda_sq_pt()=&Global_Physical_Variables::Lambda_sq;
542  }
543 
544  // Refine uniformly
545  mesh_pt()->refine_uniformly();
546 
547  // Assign equation numbers
548  assign_eqn_numbers();
549 
550 } // end constructor
551 
552 
553 
554 //=====================start_doc====================================
555 /// Doc the solution
556 //==================================================================
557 template<class ELEMENT>
559  DocInfo& doc_info)
560 {
561 
562  ofstream some_file, some_file2;
563  char filename[100];
564 
565  // Number of plot points
566  unsigned npts;
567  npts=5;
568 
569  // Output shape of deformed body
570  //------------------------------
571  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
572  doc_info.number());
573  some_file.open(filename);
574  mesh_pt()->output(some_file,npts);
575  some_file.close();
576 
577 
578  // Write trace file
579  //-----------------
580 
581  // Get position on IC object (exact solution)
582  Vector<double> r_exact(2);
583  Vector<double> xi(2);
584  xi[0]=1.0;
585  xi[1]=0.0;
586  IC_geom_object_pt->position(xi,r_exact);
587 
588  // Exact outer radius for linear elasticity
589  double exact_r=r_exact[0];
590 
591  // Add to trace file
592  Trace_file << time_pt()->time() << " "
593  << exact_r << " ";
594 
595  // Doc radii of control nodes
596  unsigned ntrace_node=Trace_node_pt.size();
597  for (unsigned j=0;j<ntrace_node;j++)
598  {
599  Trace_file << sqrt(pow(Trace_node_pt[j]->x(0),2)+
600  pow(Trace_node_pt[j]->x(1),2)) << " ";
601  }
602  Trace_file << std::endl;
603 
604 
605  // Get displacement as a function of the radial coordinate
606  //--------------------------------------------------------
607  // along boundary 0
608  //-----------------
609  {
610  // Number of elements along boundary 0:
611  unsigned nelem=mesh_pt()->nboundary_element(0);
612 
613  // Open files
614  sprintf(filename,"%s/displ_along_line%i.dat",doc_info.directory().c_str(),
615  doc_info.number());
616  some_file.open(filename);
617 
618  ofstream some_file2;
619  sprintf(filename,"%s/exact_displ_along_line%i.dat",
620  doc_info.directory().c_str(),
621  doc_info.number());
622  some_file2.open(filename);
623 
624  Vector<double> s(2);
625  Vector<double> x(2);
626  Vector<double> dxdt(2);
627  Vector<double> xi(2);
628  Vector<double> r_exact(2);
629  Vector<double> v_exact(2);
630 
631  for (unsigned e=0;e<nelem;e++)
632  {
633  some_file << "ZONE " << std::endl;
634  some_file2 << "ZONE " << std::endl;
635 
636  for (unsigned i=0;i<npts;i++)
637  {
638  // Move along bottom edge of element
639  s[0]=-1.0+2.0*double(i)/double(npts-1);
640  s[1]=-1.0;
641 
642  // Get pointer to element
643  SolidFiniteElement* el_pt=dynamic_cast<SolidFiniteElement*>
644  (mesh_pt()->boundary_element_pt(0,e));
645 
646  // Get Lagrangian coordinate
647  el_pt->interpolated_xi(s,xi);
648 
649  // Get Eulerian coordinate
650  el_pt->interpolated_x(s,x);
651 
652  // Get velocity
653  el_pt->interpolated_dxdt(s,1,dxdt);
654 
655  // Get exact Eulerian position
656  IC_geom_object_pt->position(xi,r_exact);
657 
658  // Get exact velocity
659  IC_geom_object_pt->veloc(xi,v_exact);
660 
661  // Plot radial distance and displacement
662  some_file << xi[0] << " " << x[0]-xi[0] << " "
663  << dxdt[0] << std::endl;
664 
665  some_file2 << xi[0] << " " << r_exact[0]-xi[0] << " "
666  << v_exact[0] << std::endl;
667 
668  }
669  }
670  some_file.close();
671  some_file2.close();
672 
673  } // end line output
674 
675 
676  // Get displacement (exact and computed) throughout domain
677  //--------------------------------------------------------
678  {
679  // Number of elements
680  unsigned nelem=mesh_pt()->nelement();
681 
682  // Open files
683  sprintf(filename,"%s/displ%i.dat",doc_info.directory().c_str(),
684  doc_info.number());
685  some_file.open(filename);
686 
687  sprintf(filename,"%s/exact_displ%i.dat",
688  doc_info.directory().c_str(),
689  doc_info.number());
690  some_file2.open(filename);
691 
692  Vector<double> s(2);
693  Vector<double> x(2);
694  Vector<double> dxdt(2);
695  Vector<double> xi(2);
696  Vector<double> r_exact(2);
697  Vector<double> v_exact(2);
698 
699  for (unsigned e=0;e<nelem;e++)
700  {
701  some_file << "ZONE I=" << npts << ", J="<< npts << std::endl;
702  some_file2 << "ZONE I=" << npts << ", J="<< npts << std::endl;
703 
704  for (unsigned i=0;i<npts;i++)
705  {
706  s[0]=-1.0+2.0*double(i)/double(npts-1);
707  for (unsigned j=0;j<npts;j++)
708  {
709  s[1]=-1.0+2.0*double(j)/double(npts-1);
710 
711  // Get pointer to element
712  SolidFiniteElement* el_pt=dynamic_cast<SolidFiniteElement*>
713  (mesh_pt()->element_pt(e));
714 
715  // Get Lagrangian coordinate
716  el_pt->interpolated_xi(s,xi);
717 
718  // Get Eulerian coordinate
719  el_pt->interpolated_x(s,x);
720 
721  // Get velocity
722  el_pt->interpolated_dxdt(s,1,dxdt);
723 
724  // Get exact Eulerian position
725  IC_geom_object_pt->position(xi,r_exact);
726 
727  // Get exact velocity
728  IC_geom_object_pt->veloc(xi,v_exact);
729 
730  // Plot Lagrangian position, displacement and velocity
731  some_file << xi[0] << " " << xi[1] << " "
732  << x[0]-xi[0] << " " << x[1]-xi[1] << " "
733  << dxdt[0] << " " << dxdt[1] << std::endl;
734 
735 
736  some_file2 << xi[0] << " " << xi[1] << " "
737  << r_exact[0]-xi[0] << " " << r_exact[1]-xi[1] << " "
738  << v_exact[0] << " " << v_exact[1] << std::endl;
739 
740  }
741  }
742  }
743  some_file.close();
744  some_file2.close();
745  }
746 
747 
748 }
749 
750 
751 
752 
753 //=====================start_run====================================
754 /// Run the problem: Pass number of timesteps to be performed.
755 //==================================================================
756 template<class ELEMENT>
757 void DiskOscillationProblem<ELEMENT>::run(const unsigned& nstep)
758 {
759 
760  // Output
761  DocInfo doc_info;
762 
763  // Output directory
764  doc_info.set_directory("RESLT");
765 
766  // Open trace file
767  char filename[100];
768  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
769  Trace_file.open(filename);
770 
771  // Initialise time
772  double time0=1.0;
773  time_pt()->time()=time0;
774 
775  // Set timestep
776  double dt=0.01;
777  time_pt()->initialise_dt(dt);
778 
779  // Create geometric object that specifies the initial conditions:
780 
781  // Amplitude of the oscillation
782  double ampl=0.005;
783 
784  // Build the GeomObject
785  IC_geom_object_pt=new AxisymOscillatingDisk(ampl,
787  time_stepper_pt());
788 
789  // Turn into object that specifies the initial conditions:
790  SolidInitialCondition* IC_pt = new SolidInitialCondition(IC_geom_object_pt);
791 
792  // Assign initial condition
793  SolidMesh::Solid_IC_problem.set_newmark_initial_condition_consistently(
794  this,mesh_pt(),time_stepper_pt(),IC_pt,dt,
796 
797  // Doc initial state
798  doc_solution(doc_info);
799  doc_info.number()++;
800 
801  //Timestepping loop
802  for(unsigned i=0;i<nstep;i++)
803  {
804  unsteady_newton_solve(dt);
805  doc_solution(doc_info);
806  doc_info.number()++;
807  }
808 
809 } // end of run
810 
811 
812 
813 
814 
815 
816 //========start_main====================================================
817 /// Driver for disk oscillation problem
818 //======================================================================
819 int main(int argc, char* argv[])
820 {
821 
822  // Store command line arguments
823  CommandLineArgs::setup(argc,argv);
824 
825  // If there's a command line argument run the validation (i.e. do only
826  // 10 timesteps); otherwise do a few cycles
827  unsigned nstep=1000;
828  if (CommandLineArgs::Argc!=1)
829  {
830  nstep=10;
831  }
832 
833  // Hookean constitutive equations
835  new GeneralisedHookean(&Global_Physical_Variables::Nu);
836 
837  //Set up the problem
839 
840  //Run the simulation
841  problem.run(nstep);
842 
843 } // end of main
844 
845 
846 
847 
848 
849 
850 
851 
AxisymOscillatingDisk(const double &ampl, const double &nu, TimeStepper *time_stepper_pt)
Constructor: 2 Lagrangian coordinate, 2 Eulerian coords. Pass amplitude of oscillation, Poisson ratio nu, and pointer to global timestepper.
void actions_after_newton_solve()
Update function (empty)
void accel(const Vector< double > &xi, Vector< double > &accel)
Parametrised acceleration on object at current time: accel = d^2 r(xi)/dt^2.
double Lambda_sq
Timescale ratio.
Vector< Node * > Trace_node_pt
Vector of pointers to nodes whose position we're tracing.
AxisymOscillatingDisk * IC_geom_object_pt
Geometric object that specifies the initial conditions.
void position(const Vector< double > &xi, Vector< double > &r) const
Position vector at Lagrangian coordinate xi at present time.
void run(const unsigned &nstep)
Run the problem: Pass number of timesteps to be performed.
Axisymmetrially oscillating disk with displacement field according to linear elasticity.
double multiplier(const Vector< double > &xi)
Multiplier for inertia terms (needed for consistent assignment of initial conditions in Newmark schem...
void dposition_dt(const Vector< double > &xi, const unsigned &j, Vector< double > &drdt)
Parametrised j-th time-derivative on object at current time: .
ConstitutiveLaw * Constitutive_law_pt
Pointer to constitutive law.
DiskOscillationProblem()
Constructor.
void actions_before_newton_solve()
Update function (empty)
void doc_solution(DocInfo &doc_info)
Doc the solution.
void veloc(const Vector< double > &xi, Vector< double > &veloc)
Parametrised velocity on object at current time: veloc = d r(xi)/dt.
~AxisymOscillatingDisk()
Destructor (empty)
double T
Period of oscillation.
double Nu
Poisson ratio nu.
int main(int argc, char *argv[])
Driver for disk oscillation problem.
static void residual_for_dispersion(const Vector< double > &param, const Vector< double > &omega, Vector< double > &residual)
Residual of dispersion relation for use in black-box Newton method which requires global (or static) ...
Problem class to simulate small-amplitude oscillations of a circular disk.
double Omega
Eigenfrequency.
ofstream Trace_file
Trace file.
double Nu
Poisson's ratio.
double Ampl
Amplitude of oscillation.
ElasticRefineableQuarterCircleSectorMesh< ELEMENT > * mesh_pt()
Access function for the solid mesh.