lin_unsteady_ring.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 for small amplitude ring oscillations
31 
32 //OOMPH-LIB includes
33 #include "generic.h"
34 #include "beam.h"
35 #include "meshes/one_d_lagrangian_mesh.h"
36 
37 using namespace std;
38 
39 using namespace oomph;
40 
41 ////////////////////////////////////////////////////////////////////////
42 ////////////////////////////////////////////////////////////////////////
43 ////////////////////////////////////////////////////////////////////////
44 
45 
46 //====start_of_namespace============================
47 /// Namespace for physical parameters
48 //==================================================
49 namespace Global_Physical_Variables
50 {
51 
52  /// Flag for long/short run: Default = perform long run
53  unsigned Long_run_flag=1;
54 
55  /// \short Flag for fixed timestep: Default = fixed timestep
56  unsigned Fixed_timestep_flag=1;
57 
58  /// \short Boolean flag to decide if to set IC for Newmark
59  /// directly or consistently : No Default
61 
62 } // end of namespace
63 
64 /////////////////////////////////////////////////////////////////////
65 /////////////////////////////////////////////////////////////////////
66 /////////////////////////////////////////////////////////////////////
67 
68 
69 
70 //==start_of_problem_class==============================================
71 /// Oscillating ring problem: Compare small-amplitude oscillations
72 /// against analytical solution of the linearised equations.
73 //======================================================================
74 template<class ELEMENT, class TIMESTEPPER>
75 class ElasticRingProblem : public Problem
76 {
77 
78 public:
79 
80  /// \short Constructor: Number of elements, length of domain, flag for
81  /// setting Newmark IC directly or consistently
82  ElasticRingProblem(const unsigned &N, const double &L);
83 
84  /// Access function for the mesh
85  OneDLagrangianMesh<ELEMENT>* mesh_pt()
86  {
87  return dynamic_cast<OneDLagrangianMesh<ELEMENT>*>(Problem::mesh_pt());
88  }
89 
90  /// Update function is empty
92 
93  /// Update function is empty
95 
96  /// Doc solution
97  void doc_solution(DocInfo& doc_info);
98 
99  /// Do unsteady run
100  void unsteady_run();
101 
102 private:
103 
104  /// Length of domain (in terms of the Lagrangian coordinates)
105  double Length;
106 
107  /// \short In which element are we applying displacement control?
108  /// (here only used for doc of radius)
110 
111  /// At what local coordinate are we applying displacement control?
112  Vector<double> S_displ_control;
113 
114  /// Pointer to geometric object that represents the undeformed shape
115  GeomObject* Undef_geom_pt;
116 
117  /// \short Pointer to object that specifies the initial condition
118  SolidInitialCondition* IC_pt;
119 
120  /// Trace file for recording control data
121  ofstream Trace_file;
122 }; // end of problem class
123 
124 
125 
126 
127 //===start_of_constructor===============================================
128 /// Constructor for elastic ring problem
129 //======================================================================
130 template<class ELEMENT,class TIMESTEPPER>
132 (const unsigned& N, const double& L)
133  : Length(L)
134 {
135 
136  //Allocate the timestepper -- This constructs the time object as well
137  add_time_stepper_pt(new TIMESTEPPER());
138 
139  // Undeformed beam is an elliptical ring
140  Undef_geom_pt=new Ellipse(1.0,1.0);
141 
142  //Now create the (Lagrangian!) mesh
143  Problem::mesh_pt() = new OneDLagrangianMesh<ELEMENT>(
144  N,L,Undef_geom_pt,Problem::time_stepper_pt());
145 
146  // Boundary condition:
147 
148  // Bottom:
149  unsigned ibound=0;
150  // No vertical displacement
151  mesh_pt()->boundary_node_pt(ibound,0)->pin_position(1);
152  // Zero slope: Pin type 1 dof for displacement direction 0
153  mesh_pt()->boundary_node_pt(ibound,0)->pin_position(1,0);
154 
155  // Top:
156  ibound=1;
157  // No horizontal displacement
158  mesh_pt()->boundary_node_pt(ibound,0)->pin_position(0);
159  // Zero slope: Pin type 1 dof for displacement direction 1
160  mesh_pt()->boundary_node_pt(ibound,0)->pin_position(1,1);
161 
162 
163  // Resize vector of local coordinates for control displacement
164  // (here only used to identify the point whose displacement we're
165  // tracing)
166  S_displ_control.resize(1);
167 
168  // Complete build of all elements so they are fully functional
169  // -----------------------------------------------------------
170 
171  // Find number of elements in mesh
172  unsigned Nelement = mesh_pt()->nelement();
173 
174  // Loop over the elements to set pointer to undeformed wall shape
175  for(unsigned i=0;i<Nelement;i++)
176  {
177  // Cast to proper element type
178  ELEMENT *elem_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
179 
180  // Assign the undeformed surface
181  elem_pt->undeformed_beam_pt() = Undef_geom_pt;
182  }
183 
184  // Establish control displacment: (even though no displacement
185  // control is applied we still want to doc the displacement at the same point)
186 
187  // Choose element: (This is the last one)
188  Displ_control_elem_pt=dynamic_cast<ELEMENT*>(
189  mesh_pt()->element_pt(Nelement-1));
190 
191  // Fix/doc the displacement in the vertical (1) direction at right end of
192  // the control element
193  S_displ_control[0]=1.0;
194 
195  // Do equation numbering
196  cout << "# of dofs " << assign_eqn_numbers() << std::endl;
197 
198  // Geometric object that specifies the initial conditions
199  double eps_buckl=1.0e-2;
200  double HoR=dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(0))->h();
201  unsigned n_buckl=2;
202  unsigned imode=2;
203  GeomObject* ic_geom_object_pt=
204  new PseudoBucklingRing(eps_buckl,HoR,n_buckl,imode,
205  Problem::time_stepper_pt());
206 
207  // Setup object that specifies the initial conditions:
208  IC_pt = new SolidInitialCondition(ic_geom_object_pt);
209 
210 } // end of constructor
211 
212 
213 //===start_of_doc_solution================================================
214 /// Document solution
215 //========================================================================
216 template<class ELEMENT, class TIMESTEPPER>
218  DocInfo& doc_info)
219 {
220 
221  cout << "Doc-ing step " << doc_info.number()
222  << " for time " << time_stepper_pt()->time_pt()->time() << std::endl;
223 
224 
225  // Loop over all elements to get global kinetic and potential energy
226  unsigned Nelem=mesh_pt()->nelement();
227  double global_kin=0;
228  double global_pot=0;
229  double pot,kin;
230  for (unsigned ielem=0;ielem<Nelem;ielem++)
231  {
232  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(ielem))->get_energy(pot,kin);
233  global_kin+=kin;
234  global_pot+=pot;
235  }
236 
237 
238  // Control displacement for initial condition object
239  Vector<double> xi_ctrl(1);
240  Vector<double> posn_ctrl(2);
241 
242  // Lagrangian coordinate of control point
243  xi_ctrl[0]=Displ_control_elem_pt->interpolated_xi(S_displ_control,0);
244 
245  // Get position
246  IC_pt->geom_object_pt()->position(xi_ctrl,posn_ctrl);
247 
248  // Write trace file: Time, control position, energies
249  Trace_file << time_pt()->time() << " "
250  << Displ_control_elem_pt->interpolated_x(S_displ_control,1)
251  << " " << global_pot << " " << global_kin
252  << " " << global_pot + global_kin
253  << " " << posn_ctrl[1]
254  << std::endl;
255 
256 
257  ofstream some_file;
258  char filename[100];
259 
260  // Number of plot points
261  unsigned npts=5;
262 
263  // Output solution
264  sprintf(filename,"%s/ring%i.dat",doc_info.directory().c_str(),
265  doc_info.number());
266  some_file.open(filename);
267  mesh_pt()->output(some_file,npts);
268  some_file.close();
269 
270  // Loop over all elements do dump out previous solutions
271  unsigned nsteps=time_stepper_pt()->nprev_values();
272  for (unsigned t=0;t<=nsteps;t++)
273  {
274  sprintf(filename,"%s/ring%i-%i.dat",doc_info.directory().c_str(),
275  doc_info.number(),t);
276  some_file.open(filename);
277  unsigned Nelem=mesh_pt()->nelement();
278  for (unsigned ielem=0;ielem<Nelem;ielem++)
279  {
280  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(ielem))->
281  output(t,some_file,npts);
282  }
283  some_file.close();
284  }
285 
286  // Output for initial condition object
287  sprintf(filename,"%s/ic_ring%i.dat",doc_info.directory().c_str(),
288  doc_info.number());
289  some_file.open(filename);
290 
291  unsigned nplot=1+(npts-1)*mesh_pt()->nelement();
292  Vector<double> xi(1);
293  Vector<double> posn(2);
294  Vector<double> veloc(2);
295  Vector<double> accel(2);
296 
297  for (unsigned iplot=0;iplot<nplot;iplot++)
298  {
299  xi[0]=Length/double(nplot-1)*double(iplot);
300 
301  IC_pt->geom_object_pt()->position(xi,posn);
302  IC_pt->geom_object_pt()->dposition_dt(xi,1,veloc);
303  IC_pt->geom_object_pt()->dposition_dt(xi,2,accel);
304 
305  some_file << posn[0] << " " << posn[1] << " "
306  << xi[0] << " "
307  << veloc[0] << " " << veloc[1] << " "
308  << accel[0] << " " << accel[1] << " "
309  << sqrt(pow(posn[0],2)+pow(posn[1],2)) << " "
310  << sqrt(pow(veloc[0],2)+pow(veloc[1],2)) << " "
311  << sqrt(pow(accel[0],2)+pow(accel[1],2)) << " "
312  << std::endl;
313  }
314 
315  some_file.close();
316 } // end of doc solution
317 
318 
319 
320 //===start_of_unsteady_run=================================================
321 /// Solver loop to perform unsteady run
322 //=========================================================================
323 template<class ELEMENT,class TIMESTEPPER>
325 {
326 
327  /// Label for output
328  DocInfo doc_info;
329 
330  // Output directory
331  doc_info.set_directory("RESLT");
332 
333  // Step number
334  doc_info.number()=0;
335 
336  // Set up trace file
337  char filename[100];
338  sprintf(filename,"%s/trace_ring.dat",doc_info.directory().c_str());
339  Trace_file.open(filename);
340 
341  Trace_file << "VARIABLES=\"time\",\"R<sub>ctrl</sub>\",\"E<sub>pot</sub>\"";
342  Trace_file << ",\"E<sub>kin</sub>\",\"E<sub>kin</sub>+E<sub>pot</sub>\"";
343  Trace_file << ",\"<sub>exact</sub>R<sub>ctrl</sub>\""
344  << std::endl;
345 
346  // Number of steps
347  unsigned nstep=600;
348  if (Global_Physical_Variables::Long_run_flag==0) {nstep=10;}
349 
350  // Initial timestep
351  double dt=1.0;
352 
353  // Ratio for timestep reduction
354  double timestep_ratio=1.0;
355  if (Global_Physical_Variables::Fixed_timestep_flag==0) {timestep_ratio=0.995;}
356 
357  // Number of previous timesteps stored
358  unsigned ndt=time_stepper_pt()->time_pt()->ndt();
359 
360  // Setup vector of "previous" timesteps
361  Vector<double> dt_prev(ndt);
362  dt_prev[0]=dt;
363  for (unsigned i=1;i<ndt;i++)
364  {
365  dt_prev[i]=dt_prev[i-1]/timestep_ratio;
366  }
367 
368  // Initialise the history of previous timesteps
369  time_pt()->initialise_dt(dt_prev);
370 
371  // Initialise time
372  double time0=10.0;
373  time_pt()->time()=time0;
374 
375  // Setup analytical initial condition?
377  {
378  // Note: Time has been scaled on intrinsic timescale so
379  // we don't need to specify a multiplier for the inertia
380  // terms (the default assignment of 1.0 is OK)
381  SolidMesh::Solid_IC_problem.
382  set_newmark_initial_condition_consistently(
383  this,mesh_pt(),static_cast<TIMESTEPPER*>(time_stepper_pt()),IC_pt,dt);
384  }
385  else
386  {
387  SolidMesh::Solid_IC_problem.
388  set_newmark_initial_condition_directly(
389  this,mesh_pt(),static_cast<TIMESTEPPER*>(time_stepper_pt()),IC_pt,dt);
390  }
391 
392  //Output initial data
393  doc_solution(doc_info);
394 
395  // Time integration loop
396  for(unsigned i=1;i<=nstep;i++)
397  {
398  // Solve
399  unsteady_newton_solve(dt);
400 
401  // Doc solution
402  doc_info.number()++;
403  doc_solution(doc_info);
404 
405  // Reduce timestep
406  if (time_pt()->time()<100.0) {dt=timestep_ratio*dt;}
407  }
408 
409 } // end of unsteady run
410 
411 
412 
413 //===start_of_main=====================================================
414 /// Driver for ring that performs small-amplitude oscillations
415 //=====================================================================
416 int main(int argc, char* argv[])
417 {
418 
419  // Store command line arguments
420  CommandLineArgs::setup(argc,argv);
421 
422  /// Convert command line arguments (if any) into flags:
423  if (argc==2)
424  {
425  // Nontrivial command line input: Setup Newmark IC directly
426  // (rather than consistently with PVD)
427  if (atoi(argv[1])==1)
428  {
430  cout << "Setting Newmark IC consistently" << std::endl;
431  }
432  else
433  {
435  cout << "Setting Newmark IC directly" << std::endl;
436  }
437 
438  cout << "Not enough command line arguments specified -- using defaults."
439  << std::endl;
440  } // end of 1 argument
441  else if (argc==4)
442  {
443  cout << "Three command line arguments specified:" << std::endl;
444  // Nontrivial command line input: Setup Newmark IC directly
445  // (rather than consistently with PVD)
446  if (atoi(argv[1])==1)
447  {
449  cout << "Setting Newmark IC consistently" << std::endl;
450  }
451  else
452  {
454  cout << "Setting Newmark IC directly" << std::endl;
455  }
456  // Flag for long run
458  // Flag for fixed timestep
460  } // end of 3 arguments
461  else
462  {
463  std::string error_message =
464  "Wrong number of command line arguments. Specify one or three.\n";
465  error_message += "Arg1: Long_run_flag [0/1]\n";
466  error_message += "Arg2: Impulsive_start_flag [0/1]\n";
467  error_message += "Arg3: Restart_flag [restart_file] (optional)\n";
468 
469  throw OomphLibError(error_message,
470  OOMPH_CURRENT_FUNCTION,
471  OOMPH_EXCEPTION_LOCATION);
472  } // too many arguments
473  cout << "Setting Newmark IC consistently: "
475  cout << "Long run flag: "
477  cout << "Fixed timestep flag: "
479 
480  //Length of domain
481  double L = MathematicalConstants::Pi/2.0;
482 
483  // Number of elements
484  unsigned nelem = 13;
485 
486  //Set up the problem
488  problem(nelem,L);
489 
490  // Do unsteady run
491  problem.unsteady_run();
492 
493 } // end of main
494 
495 
496 
497 
498 
499 
500 
501 
SolidInitialCondition * IC_pt
Pointer to object that specifies the initial condition.
double Length
Length of domain (in terms of the Lagrangian coordinates)
unsigned Fixed_timestep_flag
Flag for fixed timestep: Default = fixed timestep.
void unsteady_run()
Do unsteady run.
int main(int argc, char *argv[])
Driver for ring that performs small-amplitude oscillations.
ELEMENT * Displ_control_elem_pt
In which element are we applying displacement control? (here only used for doc of radius) ...
GeomObject * Undef_geom_pt
Pointer to geometric object that represents the undeformed shape.
void doc_solution(DocInfo &doc_info)
Doc solution.
OneDLagrangianMesh< ELEMENT > * mesh_pt()
Access function for the mesh.
void actions_before_newton_solve()
Update function is empty.
Vector< double > S_displ_control
At what local coordinate are we applying displacement control?
void actions_after_newton_solve()
Update function is empty.
unsigned Long_run_flag
Flag for long/short run: Default = perform long run.
ofstream Trace_file
Trace file for recording control data.
bool Consistent_newmark_ic
Boolean flag to decide if to set IC for Newmark directly or consistently : No Default.
ElasticRingProblem(const unsigned &N, const double &L)
Constructor: Number of elements, length of domain, flag for setting Newmark IC directly or consistent...