pseudo_buckling_ring.h
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 #ifndef OOMPH_PSEUDO_BUCKLING_RING_HEADER
31 #define OOMPH_PSEUDO_BUCKLING_RING_HEADER
32 
33 
34 // Config header generated by autoconfig
35 #ifdef HAVE_CONFIG_H
36  #include <oomph-lib-config.h>
37 #endif
38 
39 //oomph-lib headers
40 #include "elements.h"
41 #include "nodes.h"
42 #include "quadtree.h"
43 #include "mesh.h"
44 #include "timesteppers.h"
45 #include "geom_objects.h"
46 
47 namespace oomph
48 {
49 
50 //=========================================================================
51 /// \short Pseudo buckling ring: Circular ring deformed by the
52 /// N-th buckling mode of a thin-wall elastic ring.
53 /// \f[
54 /// x = R_0 \cos(\zeta) +
55 /// \epsilon \left( \cos(N \zeta) \cos(\zeta) - A \sin(N \zeta) \sin(\zeta)
56 /// \right) sin(2 \pi t/T)
57 /// \f]
58 /// \f[
59 /// y = R_0 \sin(\zeta) +
60 /// \epsilon \left( \cos(N \zeta) \sin(\zeta) + A \sin(N \zeta) \cos(\zeta)
61 /// \right) sin(2 \pi t/T)
62 /// \f]
63 /// where A is the ratio of the aziumuthal to the radial buckling
64 /// amplitude (A=-1/N for statically buckling rings) and epsilon
65 /// is the buckling amplitude.
66 ///
67 //=========================================================================
69 {
70 
71 public:
72 
73 
74  /// Default constructor (empty and broken)
76  {
77  throw OomphLibError(
78  "Don't call empty constructor for PseudoBucklingRing!",
79  OOMPH_CURRENT_FUNCTION,
80  OOMPH_EXCEPTION_LOCATION);
81  }
82 
83  /// \short Constructor: 1 Lagrangian coordinate, 2 Eulerian coords. Pass
84  /// buckling amplitude, ratio of of buckling amplitudes, buckling
85  /// wavenumber (as a double), undeformed ring radius (all as Data)
86  /// and pointer to global timestepper.
87  /// \code
88  /// Geom_data_pt[0]->value(0) = Eps_buckl;
89  /// Geom_data_pt[0]->value(1) = Ampl_ratio;
90  /// Geom_data_pt[0]->value(2) = Double_N_buckl;
91  /// Geom_data_pt[0]->value(3) = R_0;
92  /// Geom_data_pt[0]->value(4) = T;
93  /// \endcode
96  GeomObject(1,2,time_stepper_pt)
97  {
98 #ifdef PARANOID
99  if (geom_data_pt.size()!=1)
100  {
101  std::ostringstream error_message;
102  error_message
103  << "geom_data_pt should be of size 1, but is of size"
104  << geom_data_pt.size() << std::endl;
105 
106  throw OomphLibError(error_message.str(),
107  OOMPH_CURRENT_FUNCTION,
108  OOMPH_EXCEPTION_LOCATION);
109  }
110  if (geom_data_pt[0]->nvalue()!=5)
111  {
112  std::ostringstream error_message;
113  error_message
114  << "geom_data_pt[0] should have 5 values, but has"
115  << geom_data_pt[0]->nvalue() << std::endl;
116 
117  throw OomphLibError(error_message.str(),
118  OOMPH_CURRENT_FUNCTION,
119  OOMPH_EXCEPTION_LOCATION);
120  }
121 #endif
122  Geom_data_pt.resize(1);
123  Geom_data_pt[0]=geom_data_pt[0];
124 
125  // Data has been created externally: Must not clean up
126  Must_clean_up=false;
127  }
128 
129 
130  /// \short Constructor: 1 Lagrangian coordinate, 2 Eulerian coords. Pass
131  /// buckling amplitude, ratio of of buckling amplitudes, buckling
132  /// wavenumber, undeformed ring radius, period of osc and pointer
133  /// to global timestepper. All geometric data is pinned by default.
135  const double& ampl_ratio, const unsigned n_buckl,
136  const double& r_0,
137  const double& T,
139  GeomObject(1,2,time_stepper_pt)
140  {
141  // Number of previous timesteps that need to be stored
142  unsigned n_time=time_stepper_pt->nprev_values();
143 
144  // Setup geometric data for element: Five items of data live in
145  // in the one and only geometric data. Also store time history
146  Geom_data_pt.resize(1);
147  Geom_data_pt[0]=new Data(time_stepper_pt,5);
148 
149  // I've created the data, I need to clean up
150  Must_clean_up=true;
151 
152  for (unsigned itime=0;itime<=n_time;itime++)
153  {
154  // Buckling amplitude
155  Geom_data_pt[0]->set_value(itime,0,eps_buckl);
156  Geom_data_pt[0]->pin(0);
157 
158  // Ratio of buckling amplitudes
159  Geom_data_pt[0]->set_value(itime,1,ampl_ratio);
160  Geom_data_pt[0]->pin(1);
161 
162  // Buckling wavenumber (as double)
163  Geom_data_pt[0]->set_value(itime,2,double(n_buckl));
164  Geom_data_pt[0]->pin(2);
165 
166  // Radius of undeformed ring
167  Geom_data_pt[0]->set_value(itime,3,r_0);
168  Geom_data_pt[0]->pin(3);
169 
170  // Period of oscillation
171  Geom_data_pt[0]->set_value(itime,4,T);
172  Geom_data_pt[0]->pin(4);
173  }
174  }
175 
176 
177  /// \short Constructor: 1 Lagrangian coordinate, 2 Eulerian coords. Pass
178  /// buckling amplitude, h/R, buckling wavenumbe and pointer
179  /// to global timestepper. Other parameters get set up to represent
180  /// oscillating ring with mode imode (1 or 2). All geometric data is
181  /// pinned by default.
183  const double& HoR,
184  const unsigned& n_buckl,
185  const unsigned& imode,
187  GeomObject(1,2,time_stepper_pt)
188  {
189  // Constants in Soedel solution:
190  double K1=(pow(double(n_buckl),2)+1.0)*
191  (1.0/12.0*pow(double(n_buckl),2)*pow(HoR,2)+1.0);
192 
193  double K2oK1sq=1.0/12.0*pow(HoR,2)*pow(double(n_buckl),2)*
194  pow( (pow(double(n_buckl),2)-1.0) ,2)/
195  ( pow((pow(double(n_buckl),2)+1.0),2) *
196  pow((1.0/12.0*pow(double(n_buckl),2)*pow(HoR,2)+1.0),2) );
197 
198  // The two fundamental frequencies:
199  double omega1=sqrt(0.5*K1*(1.0+sqrt(1.0-4*K2oK1sq)));
200  double omega2=sqrt(0.5*K1*(1.0-sqrt(1.0-4*K2oK1sq)));
201 
202  // The two amplitude ratios:
203  double ampl_ratio1=
204  (double(n_buckl)*(1.0/12.0*pow(double(n_buckl),2)*pow(HoR,2)+1.0))/
205  (pow(omega1,2)-pow(double(n_buckl),2)*(1.0/12.0*pow(HoR,2)+1.0));
206 
207  double ampl_ratio2=
208  double(n_buckl)*(1.0/12.0*pow(double(n_buckl),2)*pow(HoR,2)+1.0)/
209  (pow(omega2,2)-pow(double(n_buckl),2)*(1.0/12.0*pow(HoR,2)+1.0));
210 
211  double T;
212  double ampl_ratio;
213 
214  if (n_buckl>1)
215  {
216  T=2.0*MathematicalConstants::Pi/omega2;
217  ampl_ratio=ampl_ratio2;
218  if (imode==1)
219  {
220  T= 2.0*MathematicalConstants::Pi/omega1;
221  ampl_ratio=ampl_ratio1;
222  }
223  else if (imode==2)
224  {
225  T= 2.0*MathematicalConstants::Pi/omega2;
226  ampl_ratio=ampl_ratio2;
227  }
228  else
229  {
230  oomph_info << "wrong imode " << imode << std::endl;
231  }
232  }
233  else
234  {
235  T=2.0*MathematicalConstants::Pi/omega1;
236  ampl_ratio=ampl_ratio1;
237  }
238 
239  // Number of previous timesteps that need to be stored
240  unsigned n_time=time_stepper_pt->nprev_values();
241 
242  // Setup geometric data for element: Five items of data live in
243  // in the one and only geometric data. Also store time history
244  Geom_data_pt.resize(1);
245  Geom_data_pt[0]=new Data(time_stepper_pt,5);
246 
247  // I've created the data, I need to clean up
248  Must_clean_up=true;
249 
250  for (unsigned itime=0;itime<=n_time;itime++)
251  {
252 
253  // Buckling amplitude
254  Geom_data_pt[0]->set_value(itime,0,eps_buckl);
255  Geom_data_pt[0]->pin(0);
256 
257  // Ratio of buckling amplitudes
258  Geom_data_pt[0]->set_value(itime,1,ampl_ratio);
259  Geom_data_pt[0]->pin(1);
260 
261  // Buckling wavenumber (as double)
262  Geom_data_pt[0]->set_value(itime,2,double(n_buckl));
263  Geom_data_pt[0]->pin(2);
264 
265  // Radius of undeformed ring
266  Geom_data_pt[0]->set_value(itime,3,1.0);
267  Geom_data_pt[0]->pin(3);
268 
269  // Period of oscillation
270  Geom_data_pt[0]->set_value(itime,4,T);
271  Geom_data_pt[0]->pin(4);
272  }
273  }
274 
275  /// Broken copy constructor
277  {
278  BrokenCopy::broken_copy("PseudoBucklingRing");
279  }
280 
281  /// Broken assignment operator
283  {
284  BrokenCopy::broken_assign("PseudoBucklingRing");
285  }
286 
287 
288  /// Destructor: Clean up if necessary
290  {
291  // Do I need to clean up?
292  if (Must_clean_up)
293  {
294  delete Geom_data_pt[0];
295  Geom_data_pt[0]=0;
296  }
297  }
298 
299 
300  /// Access function for buckling amplitude
301  double eps_buckl(){return Geom_data_pt[0]->value(0);}
302 
303  /// Access function for amplitude ratio
304  double ampl_ratio(){return Geom_data_pt[0]->value(1);}
305 
306  /// Access function for undeformed radius
307  double r_0(){return Geom_data_pt[0]->value(3);}
308 
309  /// Access function for period of oscillation
310  double T(){return Geom_data_pt[0]->value(4);}
311 
312  /// Access function for buckling wavenumber (as float)
313  double n_buckl_float(){return Geom_data_pt[0]->value(2);}
314 
315  /// Set buckling amplitude
316  void set_eps_buckl(const double& eps_buckl)
317  {Geom_data_pt[0]->set_value(0,eps_buckl);}
318 
319  /// \short Set amplitude ratio between radial and azimuthal
320  /// buckling displacements
321  void set_ampl_ratio(const double& ampl_ratio)
322  {Geom_data_pt[0]->set_value(1,ampl_ratio);}
323 
324  /// Set buckling wavenumber
325  void set_n_buckl(const unsigned& n_buckl)
326  {Geom_data_pt[0]->set_value(2,double(n_buckl));}
327 
328  /// Set undeformed radius of ring
329  void set_R_0(const double& r_0)
330  {Geom_data_pt[0]->set_value(3,r_0);}
331 
332  /// Set period of oscillation
333  void set_T(const double& T)
334  {Geom_data_pt[0]->set_value(4,T);}
335 
336 
337  /// \short Position Vector at Lagrangian coordinate zeta at present
338  /// time
339  void position(const Vector<double>& zeta, Vector<double>& r) const
340  {
341 #ifdef PARANOID
342  if (r.size()!=Ndim)
343  {
344  throw OomphLibError("The position vector r has the wrong dimension",
345  OOMPH_CURRENT_FUNCTION,
346  OOMPH_EXCEPTION_LOCATION);
347  }
348 #endif
349 
350  // Parameter values at present time
351  double time=Geom_object_time_stepper_pt->time_pt()->time();
352  double Eps_buckl = Geom_data_pt[0]->value(0);
353  double Ampl_ratio = Geom_data_pt[0]->value(1);
354  double double_N_buckl = Geom_data_pt[0]->value(2);
355  double R_0 = Geom_data_pt[0]->value(3);
356  double T = Geom_data_pt[0]->value(4);
357 
358 
359  // Position Vector
360  r[0]=R_0*cos(zeta[0]) +
361  Eps_buckl*( cos(double_N_buckl*zeta[0])*cos(zeta[0])-
362  Ampl_ratio*sin(double_N_buckl*zeta[0])*sin(zeta[0]) )
363  *sin(2.0*MathematicalConstants::Pi*time/T);
364 
365  r[1]=R_0*sin(zeta[0]) +
366  Eps_buckl*( cos(double_N_buckl*zeta[0])*sin(zeta[0])+
367  Ampl_ratio*sin(double_N_buckl*zeta[0])*cos(zeta[0]) )
368  *sin(2.0*MathematicalConstants::Pi*time/T);
369  }
370 
371 
372 
373 
374  ///\short Parametrised velocity on object at current time: veloc = d r(zeta)/dt.
375  void veloc(const Vector<double>& zeta, Vector<double>& veloc) // const
376  {
377 #ifdef PARANOID
378  if (veloc.size()!=Ndim)
379  {
380  throw OomphLibError(
381  "The vector veloc has the wrong size",
382  OOMPH_CURRENT_FUNCTION,
383  OOMPH_EXCEPTION_LOCATION);
384  }
385 #endif
386 
387  // Parameter values at present time
388  double time=Geom_object_time_stepper_pt->time_pt()->time();
389  double Eps_buckl = Geom_data_pt[0]->value(0);
390  double Ampl_ratio = Geom_data_pt[0]->value(1);
391  double double_N_buckl = Geom_data_pt[0]->value(2);
392  //Unused double R_0 = Geom_data_pt[0]->value(3);
393  double T = Geom_data_pt[0]->value(4);
394 
395  // Veloc
396  veloc[0]=
397  Eps_buckl*( cos(double_N_buckl*zeta[0])*cos(zeta[0])-
398  Ampl_ratio*sin(double_N_buckl*zeta[0])*sin(zeta[0]) )
400  veloc[1]=
401  Eps_buckl*( cos(double_N_buckl*zeta[0])*sin(zeta[0])+
402  Ampl_ratio*sin(double_N_buckl*zeta[0])*cos(zeta[0]) )
404  }
405 
406 
407  /// \short Parametrised acceleration on object at current time:
408  /// accel = d^2 r(zeta)/dt^2.
409  void accel(const Vector<double>& zeta, Vector<double>& accel) //const
410  {
411 #ifdef PARANOID
412  if (accel.size()!=Ndim)
413  {
414  throw OomphLibError("The vector accel has the wrong dimension",
415  OOMPH_CURRENT_FUNCTION,
416  OOMPH_EXCEPTION_LOCATION);
417  }
418 #endif
419 
420  // Parameter values at present time
421  double time=Geom_object_time_stepper_pt->time_pt()->time();
422  double Eps_buckl = Geom_data_pt[0]->value(0);
423  double Ampl_ratio = Geom_data_pt[0]->value(1);
424  double double_N_buckl = Geom_data_pt[0]->value(2);
425  //Unused double R_0 = Geom_data_pt[0]->value(3);
426  double T = Geom_data_pt[0]->value(4);
427 
428  // Accel
429  accel[0]=
430  -Eps_buckl*( cos(double_N_buckl*zeta[0])*cos(zeta[0])-
431  Ampl_ratio*sin(double_N_buckl*zeta[0])*sin(zeta[0]) )
433 
434  accel[1]=
435  -Eps_buckl*( cos(double_N_buckl*zeta[0])*sin(zeta[0])+
436  Ampl_ratio*sin(double_N_buckl*zeta[0])*cos(zeta[0]) )
438  }
439 
440 
441  /// \short Position Vector at Lagrangian coordinate zeta at discrete
442  /// previous time (t=0: present time; t>0: previous time)
443  void position(const unsigned& t, const Vector<double>& zeta,
444  Vector<double>& r) const
445  {
446 #ifdef PARANOID
447  if (r.size()!=Ndim)
448  {
449  throw OomphLibError("The position vector r has the wrong dimension",
450  OOMPH_CURRENT_FUNCTION,
451  OOMPH_EXCEPTION_LOCATION);
452  }
454  {
455  throw OomphLibError(
456  "The time value t is greater than the number of previous steps",
457  OOMPH_CURRENT_FUNCTION,
458  OOMPH_EXCEPTION_LOCATION);
459  }
460 #endif
461 
462  // Parameter values at previous time
463  double Eps_buckl = Geom_data_pt[0]->value(t,0);
464  double Ampl_ratio = Geom_data_pt[0]->value(t,1);
465  double double_N_buckl = Geom_data_pt[0]->value(t,2);
466  double R_0 = Geom_data_pt[0]->value(t,3);
467  double T = Geom_data_pt[0]->value(4);
468 
469  // Present time
470  double time=Geom_object_time_stepper_pt->time_pt()->time();
471 
472  // Recover time at previous timestep
473  for (unsigned i=0;i<t;i++)
474  {
476  }
477 
478  // Position Vector
479  r[0]=R_0*cos(zeta[0]) +
480  Eps_buckl*( cos(double_N_buckl*zeta[0])*cos(zeta[0])-
481  Ampl_ratio*sin(double_N_buckl*zeta[0])*sin(zeta[0]) )
482  *sin(2.0*MathematicalConstants::Pi*time/T);
483 
484  r[1]=R_0*sin(zeta[0]) +
485  Eps_buckl*( cos(double_N_buckl*zeta[0])*sin(zeta[0])+
486  Ampl_ratio*sin(double_N_buckl*zeta[0])*cos(zeta[0]) )
487  *sin(2.0*MathematicalConstants::Pi*time/T);
488 
489  }
490 
491 
492  /// \short j-th time-derivative on object at current time:
493  /// \f$ \frac{d^{j} r(\zeta)}{dt^j} \f$.
494  void dposition_dt(const Vector<double>& zeta, const unsigned& j,
495  Vector<double>& drdt)
496  {
497  switch (j)
498  {
499  // Current position
500  case 0:
501  position(zeta,drdt);
502  break;
503 
504  // Velocity:
505  case 1:
506  veloc(zeta,drdt);
507  break;
508 
509  // Acceleration:
510  case 2:
511  accel(zeta,drdt);
512  break;
513 
514  default:
515  std::ostringstream error_message;
516  error_message << j << "th derivative not implemented\n";
517 
518  throw OomphLibError(error_message.str(),
519  OOMPH_CURRENT_FUNCTION,
520  OOMPH_EXCEPTION_LOCATION);
521  }
522  }
523 
524 
525  /// How many items of Data does the shape of the object depend on?
526  unsigned ngeom_data() const {return Geom_data_pt.size();}
527 
528  /// \short Return pointer to the j-th Data item that the object's
529  /// shape depends on
530  Data* geom_data_pt(const unsigned& j) {return Geom_data_pt[j];}
531 
532 
533 
534 protected:
535 
536  /// \short Vector of pointers to Data items that affects the object's shape
538 
539  /// Do I need to clean up?
541 
542 };
543 
544 
545 ///////////////////////////////////////////////////////////////////////
546 ///////////////////////////////////////////////////////////////////////
547 // Pseudo buckling ring as element
548 ///////////////////////////////////////////////////////////////////////
549 ///////////////////////////////////////////////////////////////////////
550 
551 
552 //=========================================================================
553 /// \short Pseudo buckling ring: Circular ring deformed by the
554 /// N-th buckling mode of a thin-wall elastic ring.
555 /// \f[
556 /// x = R_0 \cos(\zeta) +
557 /// \epsilon \left( \cos(N \zeta) \cos(\zeta) - A \sin(N \zeta) \sin(\zeta)
558 /// \right) sin(2 \pi t/T)
559 /// \f]
560 /// \f[
561 /// y = R_0 \sin(\zeta) +
562 /// \epsilon \left( \cos(N \zeta) \sin(\zeta) + A \sin(N \zeta) \cos(\zeta)
563 /// \right) sin(2 \pi t/T)
564 /// \f]
565 /// where A is the ratio of the aziumuthal to the radial buckling
566 /// amplitude (A=-1/N for statically buckling rings) and epsilon
567 /// is the buckling amplitude.
568 /// Scale R_0 is adjusted to ensure conservation of (computational)
569 /// volume/area. This is implemented by a
570 /// pseudo-elasticity approach: The governing equation for \f$ R_0 \f$ is:
571 /// \f[
572 /// p_{ref} = R_0 - 1.0
573 /// \f]
574 /// The pointer to the reference pressure needs to be set with
575 /// reference_pressure_pt().
576 //=========================================================================
578  public PseudoBucklingRing
579 {
580  private:
581 
582  /// \short Index of the value stored in the single geometric object that has
583  /// become an unknown
585 
586  /// \short The Data object that represents the reference pressure is stored
587  /// at the location indexed by this integer in the external data storage.
589 
590  /// Pointer to the data object that represents the external reference pressure
592 
593  /// Return the local equation number of the internal geometric variable
594  inline int geometric_local_eqn()
596 
597  /// Return the local equation number of the reference pressure variable
600 
601 
602 public:
603 
604  /// \short Constructor: Build pseudo buckling ring
605  /// from doubles that describe the geometry.
607  const double& ampl_ratio, const unsigned n_buckl,
608  const double& r_0, const double& T,
610  PseudoBucklingRing(eps_buckl,ampl_ratio,n_buckl,
611  r_0,T,time_stepper_pt),
613  {
614  // Geom data for geom object has been setup (and pinned) in
615  // constructor for geometric object. Now free the scale for the half axes
616  // because we want to determine it as an unknown
617  Geom_data_pt[0]->unpin(3);
618 
619  //Record that the geometric variable is value 3 in the geometric data
621 
622  // The geometric data is internal to the element -- this
623  // ensures that any unknown pieces of geom_data get global equation
624  // numbers. There should only be one piece of internal data
625  unsigned n_geom_data=Geom_data_pt.size();
626  for (unsigned i=0;i<n_geom_data;i++) {add_internal_data(Geom_data_pt[i]);}
627  }
628 
629 
630  /// \short Constructor: Pass
631  /// buckling amplitude, h/R, buckling wavenumbe and pointer
632  /// to global timestepper. Other parameters get set up to represent
633  /// oscillating ring with mode imode (1 or 2). All geometric data is
634  /// pinned by default.
636  const double& HoR,
637  const unsigned& n_buckl,
638  const unsigned& imode,
640  PseudoBucklingRing(eps_buckl,HoR,n_buckl,imode,time_stepper_pt),
642  {
643  // Geom data for geom object has been setup (and pinned) in
644  // constructor for geometric object. Now free the scale for the half axes
645  // because we want to determine it as an unknown
646  Geom_data_pt[0]->unpin(3);
647 
648  //Record that the geometric variable is value 3 in the geometric data
650 
651  // The geometric data is internal to the element -- this
652  // ensures that any unknown pieces of geom_data get global equation
653  // numbers. There should only be one piece of internal data
654  unsigned n_geom_data=Geom_data_pt.size();
655  for (unsigned i=0;i<n_geom_data;i++) {add_internal_data(Geom_data_pt[i]);}
656  }
657 
658 
659  /// Broken copy constructor
661  {
662  BrokenCopy::broken_copy("PseudoBucklingRingElement");
663  }
664 
665  /// Broken assignment operator
667  {
668  BrokenCopy::broken_assign("PseudoBucklingRingElement");
669  }
670 
671  /// Destructor: Kill internal data and set to NULL
673  {
674  // The GeomElement's GeomData is mirrored in the element's
675  // Internal Data and therefore gets wiped in the
676  // destructor of GeneralisedElement --> No need to
677  // kill it in PseudoBucklingRing()
678  Must_clean_up=false;
679  }
680 
681 
682  /// Compute element residual Vector (wrapper)
683  inline virtual void get_residuals(Vector<double> &residuals)
684  {
685  //Create a dummy matrix
686  DenseMatrix<double> dummy(1);
687 
688  //Call the generic residuals function with flag set to 0
689  get_residuals_generic(residuals,dummy,0);
690  }
691 
692 
693  /// Compute element residual Vector and element Jacobian matrix (wrapper)
694  inline virtual void get_jacobian(Vector<double> &residuals,
695  DenseMatrix<double> &jacobian)
696  {
697  //Call the generic routine with the flag set to 1
698  get_residuals_generic(residuals,jacobian,1);
699  }
700 
701  /// \short Pointer to pressure data that is used as reference pressure
702  Data* const &reference_pressure_pt() const {return external_data_pt(0);}
703 
704  /// \short Return the reference pressure
705  double reference_pressure() const
706  {
707  //If there is no external pressure, return 0.0
708  if(External_reference_pressure_pt == 0) {return 0.0;}
709  else {return External_reference_pressure_pt->value(0);}
710  }
711 
712  /// \short Set the pressure data that is used as reference pressure
713  void set_reference_pressure_pt(Data* const &data_pt)
714  {
715  //Clear the existing external data, if there is any
717  //Set the new External reference pointer
719  //Add it to the external data
721  }
722 
723  protected:
724 
725  /// \short Compute element residual Vector (only if flag=0) and also
726  /// element Jacobian matrix (if flag=1)
727  virtual void get_residuals_generic(Vector<double> &residuals,
728  DenseMatrix<double> &jacobian,
729  unsigned flag)
730  {
731  //Initialise the residuals to zero
732  residuals.initialise(0.0);
733  //If computing the Jacobian initialise to zero
734  if(flag) {jacobian.initialise(0.0);}
735 
736  // There is only one equation, which is due to the internal degree
737  //of freedom
738  int local_eqn = geometric_local_eqn();
739 
740  //If it's not a boundary condition
741  if(local_eqn >= 0)
742  {
743  // Pseudo force balance
744  residuals[local_eqn] = reference_pressure() - (r_0()-1.0);
745 
746  // Work out jacobian: d residual[0]/d r_0
747  if(flag)
748  {
749  //The derivative wrt the internal unknown is
750  //Derivative residual w.r.t. scale
751  jacobian(local_eqn,local_eqn) = -1.0;
752 
753  int local_unknown = reference_pressure_local_eqn();
754  if(local_unknown >= 0)
755  {
756  jacobian(local_eqn,local_unknown) = 1.0;
757  }
758  }
759  }
760  }
761 
762 };
763 
764 }
765 
766 #endif
767 
768 
769 
double eps_buckl()
Access function for buckling amplitude.
A Generalised Element class.
Definition: elements.h:76
void position(const Vector< double > &zeta, Vector< double > &r) const
Position Vector at Lagrangian coordinate zeta at present time.
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data...
Definition: elements.h:268
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
virtual void get_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute element residual Vector and element Jacobian matrix (wrapper)
void set_T(const double &T)
Set period of oscillation.
void initialise(const T &val)
Initialize all values in the matrix to val.
Definition: matrices.h:523
void accel(const Vector< double > &zeta, Vector< double > &accel)
Parametrised acceleration on object at current time: accel = d^2 r(zeta)/dt^2.
virtual void get_residuals_generic(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Compute element residual Vector (only if flag=0) and also element Jacobian matrix (if flag=1) ...
void dposition_dt(const Vector< double > &zeta, const unsigned &j, Vector< double > &drdt)
j-th time-derivative on object at current time: .
void operator=(const PseudoBucklingRingElement &)
Broken assignment operator.
Data * External_reference_pressure_pt
Pointer to the data object that represents the external reference pressure.
void veloc(const Vector< double > &zeta, Vector< double > &veloc)
Parametrised velocity on object at current time: veloc = d r(zeta)/dt.
int geometric_local_eqn()
Return the local equation number of the internal geometric variable.
Pseudo buckling ring: Circular ring deformed by the N-th buckling mode of a thin-wall elastic ring...
PseudoBucklingRingElement(const double &eps_buckl, const double &HoR, const unsigned &n_buckl, const unsigned &imode, TimeStepper *time_stepper_pt)
Constructor: Pass buckling amplitude, h/R, buckling wavenumbe and pointer to global timestepper...
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object's shape depends on.
cstr elem_len * i
Definition: cfortran.h:607
~PseudoBucklingRing()
Destructor: Clean up if necessary.
PseudoBucklingRing()
Default constructor (empty and broken)
const double Pi
50 digits from maple
void operator=(const PseudoBucklingRing &)
Broken assignment operator.
char t
Definition: cfortran.h:572
Data *const & reference_pressure_pt() const
Pointer to pressure data that is used as reference pressure.
OomphInfo oomph_info
void set_reference_pressure_pt(Data *const &data_pt)
Set the pressure data that is used as reference pressure.
void set_eps_buckl(const double &eps_buckl)
Set buckling amplitude.
unsigned Internal_geometric_variable_index
Index of the value stored in the single geometric object that has become an unknown.
Vector< Data * > Geom_data_pt
Vector of pointers to Data items that affects the object's shape.
Pseudo buckling ring: Circular ring deformed by the N-th buckling mode of a thin-wall elastic ring...
TimeStepper * Geom_object_time_stepper_pt
Timestepper (used to handle access to geometry at previous timesteps)
Definition: geom_objects.h:407
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Definition: elements.cc:293
void position(const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
Position Vector at Lagrangian coordinate zeta at discrete previous time (t=0: present time; t>0: prev...
unsigned ngeom_data() const
How many items of Data does the shape of the object depend on?
void flush_external_data()
Flush all external data.
Definition: elements.cc:365
PseudoBucklingRingElement(const PseudoBucklingRingElement &dummy)
Broken copy constructor.
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
unsigned External_reference_pressure_index
The Data object that represents the reference pressure is stored at the location indexed by this inte...
virtual void get_residuals(Vector< double > &residuals)
Compute element residual Vector (wrapper)
PseudoBucklingRing(const double &eps_buckl, const double &ampl_ratio, const unsigned n_buckl, const double &r_0, const double &T, TimeStepper *time_stepper_pt)
Constructor: 1 Lagrangian coordinate, 2 Eulerian coords. Pass buckling amplitude, ratio of of bucklin...
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
double & dt(const unsigned &t=0)
Definition: timesteppers.h:137
virtual ~PseudoBucklingRingElement()
Destructor: Kill internal data and set to NULL.
double r_0()
Access function for undeformed radius.
void set_n_buckl(const unsigned &n_buckl)
Set buckling wavenumber.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
double ampl_ratio()
Access function for amplitude ratio.
void initialise(const _Tp &__value)
Iterate over all values and set to the desired value.
Definition: Vector.h:171
PseudoBucklingRingElement(const double &eps_buckl, const double &ampl_ratio, const unsigned n_buckl, const double &r_0, const double &T, TimeStepper *time_stepper_pt)
Constructor: Build pseudo buckling ring from doubles that describe the geometry.
Data *& external_data_pt(const unsigned &i)
Return a pointer to i-th external data object.
Definition: elements.h:662
bool Must_clean_up
Do I need to clean up?
double T()
Access function for period of oscillation.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
PseudoBucklingRing(const PseudoBucklingRing &node)
Broken copy constructor.
PseudoBucklingRing(const double &eps_buckl, const double &HoR, const unsigned &n_buckl, const unsigned &imode, TimeStepper *time_stepper_pt)
Constructor: 1 Lagrangian coordinate, 2 Eulerian coords. Pass buckling amplitude, h/R...
PseudoBucklingRing(const Vector< Data * > &geom_data_pt, TimeStepper *time_stepper_pt)
Constructor: 1 Lagrangian coordinate, 2 Eulerian coords. Pass buckling amplitude, ratio of of bucklin...
double reference_pressure() const
Return the reference pressure.
void set_R_0(const double &r_0)
Set undeformed radius of ring.
void set_ampl_ratio(const double &ampl_ratio)
Set amplitude ratio between radial and azimuthal buckling displacements.
int reference_pressure_local_eqn()
Return the local equation number of the reference pressure variable.
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
Definition: elements.cc:66
double n_buckl_float()
Access function for buckling wavenumber (as float)
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:291
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
unsigned Ndim
Number of Eulerian coordinates.
Definition: geom_objects.h:403
int external_local_eqn(const unsigned &i, const unsigned &j)
Return the local equation number corresponding to the j-th value stored at the i-th external data...
Definition: elements.h:313