periodic_orbit_handler.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: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #ifndef OOMPH_PERIODIC_ORBIT_HANDLER_CLASS_HEADER
31 #define OOMPH_PERIODIC_ORBIT_HANDLER_CLASS_HEADER
32 
33 // Config header generated by autoconfig
34 #ifdef HAVE_CONFIG_H
35  #include <oomph-lib-config.h>
36 #endif
37 
38 //OOMPH-LIB headers
39 #include "matrices.h"
40 #include "linear_solver.h"
42 #include "problem.h"
43 #include "assembly_handler.h"
45 #include "../meshes/one_d_mesh.template.h"
46 #include "../meshes/one_d_mesh.template.cc"
47 
48 namespace oomph
49 {
50 
51  class PeriodicOrbitEquations;
52 
53  class PeriodicOrbitAssemblyHandlerBase;
54 
55 //====================================================================
56 /// \short Timestepper used to calculate periodic orbits directly. It's
57 /// not really a "timestepper" per se, but represents the time storage
58 /// and means of calculating time-derivatives given the underlying
59 /// discretisation.
60 //====================================================================
62 {
63  friend class PeriodicOrbitEquations;
64 
65  public:
66 
67  ///Constructor for the case when we allow adaptive timestepping
68  PeriodicOrbitTimeDiscretisation(const unsigned &n_tstorage) :
69  TimeStepper(n_tstorage,1)
70  {
71  Type="PeriodicOrbitTimeDiscretisation";
72  }
73 
74 
75  /// Broken copy constructor
77  {
78  BrokenCopy::broken_copy("PeriodicOrbitTimeDiscretisation");
79  }
80 
81  /// Broken assignment operator
83  {
84  BrokenCopy::broken_assign("PeriodicOrbitTimeDiscretisation");
85  }
86 
87  ///Return the actual order of the scheme
88  unsigned order() const {return 1;}
89 
90  /// \short Broken initialisation the time-history for the Data values
91  /// corresponding to an impulsive start.
92  void assign_initial_values_impulsive(Data* const &data_pt)
93  {
94  throw OomphLibError(
95  "Cannot perform impulsive start for PeriodicOrbitTimeDiscretisation",
96  OOMPH_CURRENT_FUNCTION,
97  OOMPH_EXCEPTION_LOCATION);
98  }
99 
100  /// \short Broken initialisation of
101  /// the positions for the node corresponding to an impulsive start
103  {
104  throw OomphLibError(
105  "Cannot perform impulsive start for PeriodicOrbitTimeDiscretisation",
106  OOMPH_CURRENT_FUNCTION,
107  OOMPH_EXCEPTION_LOCATION);
108  }
109 
110 
111  /// \short Typedef for function that returns the (scalar) initial
112  /// value at a given value of the continuous time t.
113  typedef double (*InitialConditionFctPt)(const double& t);
114 
115  /// \short Initialise the time-history for the Data values,
116  /// corresponding to given time history, specified by
117  /// Vector of function pointers.
118  void assign_initial_data_values(Data* const &data_pt,
120  initial_value_fct)
121  {
122  // The time history stores the previous function values
123  unsigned n_time_value = ntstorage();
124 
125  //Find number of values stored
126  unsigned n_value = data_pt->nvalue();
127 
128  //Loop over current and stored timesteps
129  for(unsigned t=0;t<n_time_value;t++)
130  {
131 
132  // Get corresponding continous time
133  double time=Time_pt->time(t);
134 
135  //Loop over values
136  for(unsigned j=0;j<n_value;j++)
137  {
138  data_pt->set_value(t,j,initial_value_fct[j](time));
139  }
140  }
141  }
142 
143  /// Broken shifting of time values
144  void shift_time_values(Data* const &data_pt)
145  {
146  throw OomphLibError(
147  "Cannot shift time values for PeriodicOrbitTimeDiscretisation",
148  OOMPH_CURRENT_FUNCTION,
149  OOMPH_EXCEPTION_LOCATION);
150  }
151 
152  /// Broken shifting of time positions
153  void shift_time_positions(Node* const &node_pt)
154  {
155  throw OomphLibError(
156  "Cannot shift time positions for PeriodicOrbitTimeDiscretisation",
157  OOMPH_CURRENT_FUNCTION,
158  OOMPH_EXCEPTION_LOCATION);
159  }
160 
161  /// Set the weights
162  void set_weights()
163  {
164 
165  }
166 
167  /// Number of previous values available.
168  unsigned nprev_values() const{return ntstorage();}
169 
170  /// Number of timestep increments that need to be stored by the scheme
171  unsigned ndt() const {return ntstorage();}
172 
173 };
174 
175 
176 
177 
178 
179 //Special element for integrating the residuals over one period
180 class PeriodicOrbitEquations : public virtual FiniteElement
181 {
182  //Storage for the total number of time variables
183  unsigned Ntstorage;
184 
185  //Pointer to the global variable that represents the frequency
186  double *Omega_pt;
187 
188  /// Pointer to global time.
190 
191 public:
192 
193 
194  //Constructor, do nothing
196 
197  /// Broken copy constructor
199  {
200  BrokenCopy::broken_copy("PeriodicOrbitEquations");
201  }
202 
203  /// Broken assignment operator
204 //Commented out broken assignment operator because this can lead to a conflict warning
205 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
206 //realise that two separate implementations of the broken function are the same and so,
207 //quite rightly, it shouts.
208  /*void operator=(const PeriodicOrbitEquations&)
209  {
210  BrokenCopy::broken_assign("PeriodicOrbitEquations");
211  }*/
212 
213  ///Set the pointer to the frequency
214  double* &omega_pt() {return Omega_pt;}
215 
216  ///Return the frequency
217  double omega() {return *Omega_pt;}
218 
219  ///Set the total number of time storage values
220  void set_ntstorage(const unsigned &n_tstorage) {Ntstorage = n_tstorage;}
221 
222  /// Retun the pointer to the global time
223  Time* &time_pt() {return Time_pt;}
224 
225  /// Return the pointer to the global time (const version)
226  Time* const &time_pt() const {return Time_pt;}
227 
228  /// Return the global time, accessed via the time pointer
229  double time() const
230  {
231  //If no Time_pt, return 0.0
232  if(Time_pt==0) {return 0.0;}
233  else {return Time_pt->time();}
234  }
235 
236  /// Add the element's contribution to its residual vector (wrapper)
238  PeriodicOrbitAssemblyHandlerBase* const &assembly_handler_pt,
239  GeneralisedElement* const &elem_pt,
240  Vector<double> &residuals)
241  {
242  //Call the generic residuals function with flag set to 0
243  //using a dummy matrix argument
245  assembly_handler_pt,elem_pt,
247  }
248 
249  /// Add the element's contribution to its residual vector and
250  /// element Jacobian matrix (wrapper)
252  PeriodicOrbitAssemblyHandlerBase* const &assembly_handler_pt,
253  GeneralisedElement* const &elem_pt,
254  Vector<double> &residuals,
255  DenseMatrix<double> &jacobian)
256  {
257  //Call the generic routine with the flag set to 1
259  elem_pt,
260  residuals,jacobian,1);
261  }
262 
263 
264  void orbit_output(GeneralisedElement* const &elem_pt,
265  std::ostream &outfile,
266  const unsigned &n_plot);
267 
268 
269 protected:
270 
271  /// The routine that actually does all the work!
273  PeriodicOrbitAssemblyHandlerBase* const &assembly_handler_pt,
274  GeneralisedElement* const &elem_pt,
275  Vector<double> &residuals, DenseMatrix<double> &jacobian,
276  const unsigned& flag);
277 
278  /// \short Set the timestepper weights
279  void set_timestepper_weights(const Shape &psi, const DShape &dpsidt)
280  {
281  PeriodicOrbitTimeDiscretisation* cast_time_stepper_pt =
283 
284 
285  //Zero the timestepper weights
286  unsigned n_time_dof = cast_time_stepper_pt->ntstorage();
287  for(unsigned i=0;i<n_time_dof;i++)
288  {
289  cast_time_stepper_pt->Weight(0,i) = 0.0;
290  cast_time_stepper_pt->Weight(1,i) = 0.0;
291  }
292 
293  //Cache the frequency (timescale)
294  const double inverse_timescale = this->omega();
295  //Now set the weights
296  const unsigned n_node = this->nnode();
297 
298  //Global equation for the total number of time unknowns
299  //in the problem
300  int global_eqn;
301  for(unsigned l=0;l<n_node;l++)
302  {
303  global_eqn = this->eqn_number(this->nodal_local_eqn(l,0));
304  cast_time_stepper_pt->Weight(0,global_eqn) = psi(l);
305  cast_time_stepper_pt->Weight(1,global_eqn) =
306  dpsidt(l,0)*inverse_timescale;
307  }
308  }
309 
310  /// \short Shape/test functions and derivs w.r.t. to global coords at
311  /// local coord. s; return Jacobian of mapping
312  virtual double dshape_and_dtest_eulerian_orbit(const Vector<double> &s,
313  Shape &psi,
314  DShape &dpsidt, Shape &test,
315  DShape &dtestdt) const=0;
316 
317 
318  /// \short Shape/test functions and derivs w.r.t. to global coords at
319  /// integration point ipt; return Jacobian of mapping
320  virtual double dshape_and_dtest_eulerian_at_knot_orbit(const unsigned &ipt,
321  Shape &psi,
322  DShape &dpsidt,
323  Shape &test,
324  DShape &dtestdt)
325  const=0;
326 
327  /// \short Compute element residual Vector only (if flag=and/or element
328  /// Jacobian matrix
329 
330  //Output function
331 
332  /// \short Output function:
333  /// x,y,u or x,y,z,u at n_plot^DIM plot points
334 
335 
336 };
337 
338 
339 //======================================================================
340 /// QPoissonElement elements are linear/quadrilateral/brick-shaped
341 /// Poisson elements with isoparametric interpolation for the function.
342 //======================================================================
343 template <unsigned NNODE_1D>
345  public virtual QSpectralElement<1,NNODE_1D>,
346  public virtual RefineableQSpectralElement<1>,
347  public virtual PeriodicOrbitEquations,
348  public virtual ElementWithZ2ErrorEstimator
349 {
350 
351 
352  public:
353 
354 
355  ///\short Constructor: Call constructors for QElement and
356  /// Poisson equations
360  {}
361 
362  /// Broken copy constructor
365  {
366  BrokenCopy::broken_copy("SpectralPeriodicOrbitElement");
367  }
368 
369  /// Broken assignment operator
370  /*void operator=(const SpectralPeriodicOrbitElement<NNODE_1D>&)
371  {
372  BrokenCopy::broken_assign("SpectralPeriodicOrbitElement");
373  }*/
374 
375 
376  /// \short Required # of `values' (pinned or dofs)
377  /// at node n (only ever one dummy value, used for equation numbering)
378  /// This is also used to represent all *spatial* variables during a
379  /// temporal refinement, which is a bit naughty but it quick and dirty.
380  inline unsigned required_nvalue(const unsigned &n) const
381  {return 1;}
382 
383  /// \short Number of continuously interpolated values (1)
384  inline unsigned ncont_interpolated_values() const {return 1;}
385 
386  /// \short Return the dummy values
388  {
389  this->get_interpolated_values(0,s,value);
390  }
391 
392  ///\short Return the temporal dummy values
393  void get_interpolated_values(const unsigned &t, const Vector<double> &s,
394  Vector<double> &value)
395  {
396  value.resize(1);
397  value[0] = 0.0;
398 
399  const unsigned n_node = this->nnode();
400  Shape psi(n_node);
401  this->shape(s,psi);
402 
403  for(unsigned n=0;n<n_node;n++)
404  {
405  value[0] += this->nodal_value(t,n,0)*psi(n);
406  }
407  }
408 
409  //\short Order of recovery shape functions for Z2 error estimation:
410  unsigned nrecovery_order() {return NNODE_1D-1;}
411 
412  /// \short Number of flux terms for Z2 error estimation
413  /// This will be used to represent all spatial values,
414  unsigned num_Z2_flux_terms()
415  {return this->node_pt(0)->ntstorage();}
416 
417  /// \short Get the fluxes for the recovert
419  {
420  //Find out the number of nodes in the element
421  const unsigned n_node = this->nnode();
422 
423  //Get the shape functions
424  Shape psi(n_node);
425  DShape dpsidx(n_node,1);
426 
427  //Get the derivatives
428  (void)this->dshape_eulerian(s,psi,dpsidx);
429 
430  //Now assemble all the derivatives
431  const unsigned n_tstorage = this->node_pt(0)->ntstorage();
432 
433  //Zero the flux vector
434  for(unsigned t=0;t<n_tstorage;t++)
435  {
436  flux[t] = 0.0;
437  }
438 
439  //Loop over the nodes
440  for(unsigned n=0;n<n_node;n++)
441  {
442  const double dpsidx_ = dpsidx(n,0);
443  for(unsigned t=0;t<n_tstorage;t++)
444  {
445  flux[t] += this->nodal_value(t,n,0)*dpsidx_;
446  }
447  }
448  }
449 
450  //Number of vertex nodes in the element (always 2)
451  unsigned nvertex_node() const {return 2;}
452 
453  /// \short Pointer to the j-th vertex node in the element
454  Node* vertex_node_pt(const unsigned& j) const
456 
457 
458  //
459  /// \short Function to return the number of values
460 
461  /// \short Output function:
462  /// x,y,u or x,y,z,u
463  void output(std::ostream &outfile)
465 
466 
467  /// \short Output function:
468  /// x,y,u or x,y,z,u at n_plot^DIM plot points
469  void output(std::ostream &outfile, const unsigned &n_plot)
470  {PeriodicOrbitEquations::output(outfile,n_plot);}
471 
472 
473  /// \short C-style output function:
474  /// x,y,u or x,y,z,u
475  void output(FILE* file_pt)
477 
478 
479  /// \short C-style output function:
480  /// x,y,u or x,y,z,u at n_plot^DIM plot points
481  void output(FILE* file_pt, const unsigned &n_plot)
482  {PeriodicOrbitEquations::output(file_pt,n_plot);}
483 
484 
485 
486 protected:
487 
488 /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
489  inline double dshape_and_dtest_eulerian_orbit(
490  const Vector<double> &s, Shape &psi, DShape &dpsidt,
491  Shape &test, DShape &dtestdt) const;
492 
493 
494  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
495  /// integration point ipt. Return Jacobian.
496  inline double dshape_and_dtest_eulerian_at_knot_orbit(const unsigned& ipt,
497  Shape &psi,
498  DShape &dpsidt,
499  Shape &test,
500  DShape &dtestdt)
501  const;
502 
503 };
504 
505 
506 //Inline functions:
507 
508 
509 //======================================================================
510 /// Define the shape functions and test functions and derivatives
511 /// w.r.t. global coordinates and return Jacobian of mapping.
512 ///
513 /// Galerkin: Test functions = shape functions
514 //======================================================================
515 template<unsigned NNODE_1D>
518  const Vector<double> &s,
519  Shape &psi,
520  DShape &dpsidt,
521  Shape &test,
522  DShape &dtestdt) const
523 {
524  //Call the geometrical shape functions and derivatives
525  const double J = this->dshape_eulerian(s,psi,dpsidt);
526 
527  //Set the test functions equal to the shape functions
528  test = psi;
529  dtestdt= dpsidt;
530 
531  //Return the jacobian
532  return J;
533 }
534 
535 
536 
537 //======================================================================
538 /// Define the shape functions and test functions and derivatives
539 /// w.r.t. global coordinates and return Jacobian of mapping.
540 ///
541 /// Galerkin: Test functions = shape functions
542 //======================================================================
543 template<unsigned NNODE_1D>
546  const unsigned &ipt,
547  Shape &psi,
548  DShape &dpsidt,
549  Shape &test,
550  DShape &dtestdt) const
551 {
552  //Call the geometrical shape functions and derivatives
553  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidt);
554 
555  //Set the pointers of the test functions
556  test = psi;
557  dtestdt = dpsidt;
558 
559  //Return the jacobian
560  return J;
561 }
562 
563 
564 
565 
566 template<unsigned NNODE_1D>
568 
569 
570 //======================================================================
571 /// A special temporal mesh class
572 //=====================================================================
573 template<class ELEMENT>
575 {
576  //Let's have the periodic handler as a friend
577  template<unsigned NNODE_1D>
579 
580 public:
581 
582  ///Constructor, create a 1D mesh from 0 to 1 that is periodic
583  PeriodicOrbitTemporalMesh(const unsigned &n_element) :
584  OneDMesh<ELEMENT>(n_element,1.0),
585  RefineableOneDMesh<ELEMENT>(n_element,1.0)
586  {
587  //Make the mesh periodic by setting the LAST node to have the same data
588  //as the FIRST node
589  //Not necessarily a smart move for when doing Floquet analysis
590  this->boundary_node_pt(1,0)->
591  make_periodic(this->boundary_node_pt(0,0));
592  }
593 
594 
595  //Output the orbit for all elements in the mesh
596  void orbit_output(GeneralisedElement* const &elem_pt,
597  std::ostream &outfile, const unsigned &n_plot)
598  {
599  //Loop over all elements in the mesh
600  const unsigned n_element = this->nelement();
601  for(unsigned e=0;e<n_element;e++)
602  {
603  dynamic_cast<ELEMENT*>(this->element_pt(e))
604  ->orbit_output(elem_pt,outfile,n_plot);
605  }
606  }
607 
608  //Loop over all temporal elements and assemble their contributions
610  PeriodicOrbitAssemblyHandlerBase* const &assembly_handler_pt,
611  GeneralisedElement* const &elem_pt,
612  Vector<double> &residuals)
613  {
614  //Initialise the residuals to zero
615  residuals.initialise(0.0);
616  //Loop over all elements in the mesh
617  const unsigned n_element = this->nelement();
618  for(unsigned e=0;e<n_element;e++)
619  {
620  dynamic_cast<ELEMENT*>(this->element_pt(e))
621  ->fill_in_contribution_to_integrated_residuals(
622  assembly_handler_pt,elem_pt,residuals);
623  }
624  }
625 
626 
627  //Loop over all temporal elements and assemble their contributions
628  //and the jaobian
630  PeriodicOrbitAssemblyHandlerBase* const &assembly_handler_pt,
631  GeneralisedElement* const &elem_pt,
632  Vector<double> &residuals, DenseMatrix<double> &jacobian)
633  {
634  //Initialise the residuals to zero
635  residuals.initialise(0.0); jacobian.initialise(0.0);
636  //Loop over all elements in the mesh
637  const unsigned n_element = this->nelement();
638  for(unsigned e=0;e<n_element;e++)
639  {
640  dynamic_cast<ELEMENT*>(this->element_pt(e))
641  ->fill_in_contribution_to_integrated_jacobian(
642  assembly_handler_pt,elem_pt,residuals,jacobian);
643  }
644  }
645 
646 };
647 
648 ///===============================================================
649 /// Base class to avoid template complications
650 //===============================================================
652 {
653 public:
654  //Do nothing in constructor
656 
657  //Provide interface
658  virtual void get_dofs_for_element(GeneralisedElement* const elem_pt,
659  Vector<double> &dofs)=0;
660 
661 
662  virtual void get_previous_dofs_for_element(GeneralisedElement* const elem_pt,
663  Vector<double> &dofs)=0;
664 
665 
666  virtual void set_dofs_for_element(GeneralisedElement* const elem_pt,
667  Vector<double> const &dofs)=0;
668 
669 };
670 
671 //======================================================================
672 /// A class that is used to assemble and solve the augmented system
673 /// of equations associated with calculating periodic orbits directly
674 //========================================================================
675 template<unsigned NNODE_1D>
676 class PeriodicOrbitAssemblyHandler: public PeriodicOrbitAssemblyHandlerBase
677 {
678 private:
679 
680  ///Pointer to the timestepper
682 
683  ///Pointer to the problem
685 
686  ///Storage for mesh of temporal elements
689 
690  ///Storage for the mesh of temporal elements with a simple mesh pointer
692 
693  ///Storage for the previous solution
695 
696  ///\short Store number of degrees of freedom in the original problem
697  unsigned Ndof;
698 
699  ///Storage for number of elements in the period
701 
702  ///Storage for the number of unknown time values
703  unsigned N_tstorage;
704 
705  ///Storage for the frequency of the orbit (scaled by 2pi)
706  double Omega;
707 
708  public:
709 
710  ///Constructor, initialises values and constructs mesh of elements
712  const unsigned &n_element_in_period,
713  const DenseMatrix<double>
714  &initial_guess,
715  const double &omega) :
716  Problem_pt(problem_pt), N_element_in_period(n_element_in_period),
717  Omega(omega/(2.0*MathematicalConstants::Pi))
718  {
719  //Store the current number of degrees of freedom
720  Ndof = problem_pt->ndof();
721 
722  //Create the appropriate mesh of "1D" time elements depending on
723  //the specified spectral order
724  //The time domain runs from zero to one
725  Time_mesh_pt =
727  (n_element_in_period);
728 
730 
731  //Set the error estimator
732  Time_mesh_pt->spatial_error_estimator_pt() = new Z2ErrorEstimator;
733  Time_mesh_pt->max_permitted_error() = 1.0e-2;
734  Time_mesh_pt->min_permitted_error() = 1.0e-5;
735  Time_mesh_pt->max_refinement_level() = 10;
736 
737  //Now we need to number the mesh
738  //Dummy dof pointer
739  Vector<double*> dummy_dof_pt;
740  Time_mesh_pt->assign_global_eqn_numbers(dummy_dof_pt);
741  //Assign the local equation numbers
742  Time_mesh_pt->assign_local_eqn_numbers(false);
743 
744  //Find the number of temporal degrees of freedom
745  N_tstorage = dummy_dof_pt.size();
746 
747  //Now we need to do something clever to setup the time storage schemes
748  //and the initial values (don't quite know what that is yet)
749 
750  //Create the "fake" timestepper
752  //Loop over the temporal elements and set the pointers
753  for(unsigned e=0;e<n_element_in_period;e++)
754  {
757  (Time_mesh_pt->element_pt(e));
758 
759  //Set the time and the timestepper
760  el_pt->time_pt() = problem_pt->time_pt();
761  el_pt->time_stepper_pt() = Time_stepper_pt;
762 
763  //Set the number of temporal degrees of freedom
764  el_pt->set_ntstorage(N_tstorage);
765  //Set the frequency
766  el_pt->omega_pt() = &Omega;
767  }
768 
769  //We now need to do something much more drastic which is to loop over all
770  //our the data in the problem and change the timestepper, which is going
771  //to be a real pain when I start to worry about halo nodes, etc.
772 
773  //Will need to use the appropriate mesh-level functions that have
774  //not been written yet ..
775 
776  //Let's just break if there are submeshes
777  if(problem_pt->nsub_mesh() > 0)
778  {
779  throw OomphLibError(
780  "PeriodicOrbitHandler can't cope with submeshes yet",
781  OOMPH_CURRENT_FUNCTION,
782  OOMPH_EXCEPTION_LOCATION);
783  }
784 
785  //OK now we have only one mesh
786  unsigned n_node = problem_pt->mesh_pt()->nnode();
787  for(unsigned n=0;n<n_node;n++)
788  {
789  Node* const nod_pt = problem_pt->mesh_pt()->node_pt(n);
790  nod_pt->set_time_stepper(Time_stepper_pt,false);
791  //If the unknowns are pinned then copy the value to all values
792  unsigned n_value = nod_pt->nvalue();
793  for(unsigned i=0;i<n_value;i++)
794  {
795  if(nod_pt->is_pinned(i))
796  {
797  const unsigned n_tstorage = nod_pt->ntstorage();
798  const double value = nod_pt->value(i);
799  for(unsigned t=1;t<n_tstorage;t++)
800  {
801  nod_pt->set_value(t,i,value);
802  }
803  }
804  }
805  }
806 
807  unsigned n_element = problem_pt->mesh_pt()->nelement();
808  for(unsigned e=0;e<n_element;e++)
809  {
810  unsigned n_internal =
811  problem_pt->mesh_pt()->element_pt(e)->ninternal_data();
812  for(unsigned i=0;i<n_internal;i++)
813  {
814  //Cache the internal data
815  Data* const data_pt = problem_pt->mesh_pt()->element_pt(e)
816  ->internal_data_pt(i);
817  //and set the timestepper
818  data_pt->set_time_stepper(Time_stepper_pt,false);
819 
820  //If the unknowns are pinned then copy the value to all values
821  unsigned n_value = data_pt->nvalue();
822  for(unsigned j=0;j<n_value;j++)
823  {
824  if(data_pt->is_pinned(j))
825  {
826  const unsigned n_tstorage = data_pt->ntstorage();
827  const double value = data_pt->value(j);
828  for(unsigned t=1;t<n_tstorage;t++)
829  {
830  data_pt->set_value(t,j,value);
831  }
832  }
833  }
834  }
835  }
836 
837  //Need to reassign equation numbers so that the DOF pointer, points to the
838  //newly allocated storage
839  oomph_info << "Re-allocated " << problem_pt->assign_eqn_numbers()
840  << " equation numbers\n";
841 
842  //Now's let's add all the unknowns to the problem
843  problem_pt->Dof_pt.resize(Ndof*N_tstorage + 1);
844  //This is reasonably straight forward using pointer arithmetic
845  //but this does rely on knowing how the data is stored in the
846  //Nodes which is a little nasty
847  for(unsigned i=0;i<N_tstorage;i++)
848  {
849  unsigned offset = Ndof*i;
850  for(unsigned n=0;n<Ndof;n++)
851  {
852  problem_pt->Dof_pt[offset + n] = problem_pt->Dof_pt[n] + i;
853  }
854  }
855 
856  //Add the frequency of the orbit to the unknowns
857  problem_pt->Dof_pt[Ndof*N_tstorage] = &Omega;
858 
859  //Rebuild everything
860  problem_pt->Dof_distribution_pt->build(problem_pt->communicator_pt(),
861  Ndof*N_tstorage+1,true);
862 
863 
864  //Set initial condition of constant-ness plus wobble
865  for(unsigned i=0;i<N_tstorage;i++)
866  {
867  unsigned offset = Ndof*i;
868  for(unsigned n=0;n<Ndof;n++)
869  {
870  problem_pt->dof(offset + n) = initial_guess(i,n);//problem_pt->dof(n);
871  }
872  }
873 
874  //Set the initial unkowns to be the original problem
875  Previous_dofs.resize(Ndof*N_tstorage+1,0.0);
877 
878  //Now check everything is OK ... it seems to be
879  //std::cout << problem_pt->ndof() << "\n";
880  //Let's check it
881  //for(unsigned i=0;i<problem_pt->ndof();i++)
882  // {
883  // std::cout << i << " " << problem_pt->dof(i) << "\n";
884  //}
885 
886  }
887 
888  ///Update the previous dofs
890  {
891  for(unsigned n=0;n<Ndof*N_tstorage+1;n++)
892  {
893  Previous_dofs[n] = Problem_pt->dof(n);
894  }
895  }
896 
897  ///Return the number of degrees of freedom in the element elem_pt
898  unsigned ndof(GeneralisedElement* const &elem_pt)
899  {return ((elem_pt->ndof())*N_tstorage + 1);}
900 
901  /// \short Return the global equation number of the local unknown ieqn_local
902  ///in elem_pt.
903  unsigned long eqn_number(GeneralisedElement* const &elem_pt,
904  const unsigned &ieqn_local)
905  {
906  //Get unaugmented number of (spatial) dofs in element
907  unsigned raw_ndof = elem_pt->ndof();
908  //The final variable (the period) is stored at the end
909  if(ieqn_local == raw_ndof*N_tstorage) {return N_tstorage*Ndof;}
910  //Otherwise we need to do a little more work
911  else
912  {
913  //Now find out the time level
914  unsigned t = ieqn_local/raw_ndof;
915  //and the remainder (original eqn number)
916  unsigned raw_ieqn = ieqn_local%raw_ndof;
917  //hence calculate the global value
918  return t*Ndof + elem_pt->eqn_number(raw_ieqn);
919  }
920  }
921 
922  ///Return the contribution to the residuals of the element elem_pt
923  void get_residuals(GeneralisedElement* const &elem_pt,
924  Vector<double> &residuals)
925  {Time_mesh_pt->assemble_residuals(this,elem_pt,
926  residuals);}
927 
928 
929  //Provide interface
931  Vector<double> &dofs)
932  {
933  //Find the number of dofs in the element
934  const unsigned n_elem_dof = this->ndof(elem_pt);
935  dofs.resize(n_elem_dof);
936  //Now just get the dofs corresponding to the element's unknowns from the
937  //problem dof
938  for(unsigned i=0;i<n_elem_dof;i++)
939  {
940  dofs[i] = Problem_pt->dof(this->eqn_number(elem_pt,i));
941  }
942  }
943 
945  Vector<double> &dofs)
946  {
947  //Find the number of dofs in the element
948  const unsigned n_elem_dof = this->ndof(elem_pt);
949  dofs.resize(n_elem_dof);
950  //Now just get the dofs corresponding to the element's unknowns from the
951  //problem dof
952  for(unsigned i=0;i<n_elem_dof;i++)
953  {
954  dofs[i] = Previous_dofs[this->eqn_number(elem_pt,i)];
955  }
956  }
957 
958 
960  Vector<double> const &dofs)
961  {
962  //Find the number of dofs in the element
963  const unsigned n_elem_dof = this->ndof(elem_pt);
964  //Now just get the dofs corresponding to the element's unknowns from the
965  //problem dof
966  for(unsigned i=0;i<n_elem_dof;i++)
967  {
968  Problem_pt->dof(this->eqn_number(elem_pt,i)) = dofs[i];
969  }
970  }
971 
972 
973  /// \short Calculate the elemental Jacobian matrix "d equation
974  /// / d variable" for elem_pt.
975  void get_jacobian(GeneralisedElement* const &elem_pt,
976  Vector<double> &residuals,
977  DenseMatrix<double> &jacobian)
978  {
979  Time_mesh_pt->assemble_residuals_and_jacobian(this,elem_pt,
980  residuals,jacobian);
981  }
982 
983  /// \short Calculate all desired vectors and matrices
984  /// provided by the element elem_pt.
985  //void get_all_vectors_and_matrices(
986  // GeneralisedElement* const &elem_pt,
987  // Vector<Vector<double> >&vec, Vector<DenseMatrix<double> > &matrix) {}
988 
989  /// \short Return an unsigned integer to indicate whether the
990  /// handler is a bifurcation tracking handler. The default
991  /// is zero (not)
992  //virtual int bifurcation_type() const {return 0;}
993 
994  /// \short Return a pointer to the
995  /// bifurcation parameter in bifurcation tracking problems
996  //virtual double* bifurcation_parameter_pt() const;
997 
998  /// \short Return the eigenfunction(s) associated with the bifurcation that
999  /// has been detected in bifurcation tracking problems
1000  //virtual void get_eigenfunction(Vector<DoubleVector> &eigenfunction);
1001 
1002  ///Return the contribution to the residuals of the element elem_pt
1003  void orbit_output(std::ostream &outfile, const unsigned &n_plot)
1004  {
1005  const unsigned n_element = Problem_pt->mesh_pt()->nelement();
1006  for(unsigned e=0;e<n_element;e++)
1007  {
1008  Time_mesh_pt->orbit_output(Problem_pt->mesh_pt()->element_pt(e),
1009  outfile,n_plot);}
1010  }
1011 
1012 
1013  ///Tell me the times at which you want the solution
1015  {
1016  const unsigned n_node = Time_mesh_pt->nnode();
1017  t.resize(n_node);
1018  for(unsigned n=0;n<n_node;n++) {t[n] = Time_mesh_pt->node_pt(n)->x(0);}
1019  }
1020 
1021  ///Adapt the time mesh
1023  {
1024  //First job is to compute some sort of error measure
1025  //Then we can decide how to refine this is probably best handled
1026  //separately for now
1027 
1028  //The current plan is to copy all (locally held in the case of
1029  //distributed problem) spatial degrees of freedom into the dummy
1030  //storage of the time mesh
1031 
1032  //Probably should kick this down to the mesh level...
1033 
1034  //OK, let's do it, count up all values
1035  unsigned total_n_value=0;
1036 
1037  //Firstly the global data in the mesh
1038  unsigned n_global_data = Problem_pt->nglobal_data();
1039  for(unsigned i=0;i<n_global_data;i++)
1040  {
1041  total_n_value += Problem_pt->global_data_pt(i)->nvalue();
1042  }
1043 
1044  //Now the nodal data
1045  unsigned n_node = Problem_pt->mesh_pt()->nnode();
1046  for(unsigned n=0;n<n_node;n++)
1047  {
1048  total_n_value += Problem_pt->mesh_pt()->node_pt(n)->nvalue();
1049  SolidNode* solid_nod_pt =
1050  dynamic_cast<SolidNode*>(Problem_pt->mesh_pt()->node_pt(n));
1051  if(solid_nod_pt!=0)
1052  {
1053  total_n_value += solid_nod_pt->variable_position_pt()->nvalue();
1054  }
1055  }
1056 
1057  //Now just do the internal data
1058  unsigned n_space_element = Problem_pt->mesh_pt()->nelement();
1059  for(unsigned e=0;e<n_space_element;e++)
1060  {
1061  const unsigned n_internal_data =
1063  for(unsigned i=0;i<n_internal_data;i++)
1064  {
1065  total_n_value +=
1067  }
1068  }
1069 
1070  //Now in theory I know the total number of values in the problem
1071  //So I can create another Fake timestepper
1072  TimeStepper* fake_space_time_stepper_pt
1073  = new PeriodicOrbitTimeDiscretisation(total_n_value);
1074 
1075  //Now apply this time stepper to all time nodes
1076  unsigned n_time_node = Time_mesh_pt->nnode();
1077  for(unsigned t=0;t<n_time_node;t++) //Do include the periodic one
1078  {
1079  Time_mesh_pt->node_pt(t)->set_time_stepper(fake_space_time_stepper_pt,
1080  false);
1081  }
1082 
1083  //Now I "just" copy the values into the new storage
1084  unsigned count=0;
1085  for(unsigned i=0;i<n_global_data;i++)
1086  {
1087  Data* const glob_data_pt = Problem_pt->global_data_pt(i);
1088  const unsigned n_value = glob_data_pt->nvalue();
1089  for(unsigned j=0;j<n_value;j++)
1090  {
1091  for(unsigned t=0;t<N_tstorage;t++)
1092  {
1093  //Some heavy assumptions here about the time mesh, but that's OK
1094  //because I know exactly how it's laid out
1095  Time_mesh_pt->node_pt(t)->set_value(count,0,glob_data_pt->value(t,j));
1096  }
1097  ++count;
1098  }
1099  }
1100 
1101  //Now the nodal data
1102  for(unsigned n=0;n<n_node;n++)
1103  {
1104  Node* const nod_pt = Problem_pt->mesh_pt()->node_pt(n);
1105  const unsigned n_value = nod_pt->nvalue();
1106  for(unsigned i=0;i<n_value;i++)
1107  {
1108  for(unsigned t=0;t<N_tstorage;t++)
1109  {
1110  Time_mesh_pt->node_pt(t)->set_value(count,0,nod_pt->value(t,i));
1111  }
1112  ++count;
1113  }
1114 
1115  //Now deal with the solid node data
1116  SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
1117  if(solid_nod_pt!=0)
1118  {
1119  const unsigned n_solid_value =
1120  solid_nod_pt->variable_position_pt()->nvalue();
1121  for(unsigned i=0;i<n_solid_value;i++)
1122  {
1123  for(unsigned t=0;t<N_tstorage;t++)
1124  {
1125  Time_mesh_pt->node_pt(t)->set_value(
1126  count,0,
1127  solid_nod_pt->variable_position_pt()->value(t,i));
1128  }
1129  ++count;
1130  }
1131  }
1132  }
1133 
1134  //Now just do the internal data
1135  for(unsigned e=0;e<n_space_element;e++)
1136  {
1137  const unsigned n_internal =
1139  for(unsigned i=0;i<n_internal;i++)
1140  {
1141  Data* const internal_dat_pt =
1143 
1144  const unsigned n_value = internal_dat_pt->nvalue();
1145  for(unsigned j=0;j<n_value;j++)
1146  {
1147  for(unsigned t=0;t<N_tstorage;t++)
1148  {
1149  Time_mesh_pt->node_pt(t)->
1150  set_value(count,0,internal_dat_pt->value(t,j));
1151  }
1152  ++count;
1153  }
1154  }
1155  }
1156 
1157 
1158  //Think it's done but let's check
1159  /*{
1160  std::ofstream munge("data_remunge.dat");
1161  const unsigned n_time_element = Time_mesh_pt->nelement();
1162  for(unsigned e=0;e<n_time_element;e++)
1163  {
1164  const unsigned n_node = Time_mesh_pt->nnode();
1165  for(unsigned n=0;n<n_node;n++)
1166  {
1167  Node* const nod_pt = Time_mesh_pt->node_pt(n);
1168  munge << nod_pt->x(0) << " ";
1169  const unsigned n_space_storage = nod_pt->ntstorage();
1170  for(unsigned t=0;t<n_space_storage;t++)
1171  {
1172  munge << nod_pt->value(t,0) << " ";
1173  }
1174  munge << std::endl;
1175  }
1176  }
1177  munge.close();
1178  }*/
1179 
1180  //Ok get the elemental errors
1181  const unsigned n_time_element = Time_mesh_pt->nelement();
1182  Vector<double> elemental_error(n_time_element);
1183  Time_mesh_pt->spatial_error_estimator_pt()
1184  ->get_element_errors(Problem_pt->communicator_pt(),
1185  Basic_time_mesh_pt,elemental_error);
1186 
1187  //Let's dump it
1188  for(unsigned e=0;e<n_time_element;e++)
1189  {
1190  oomph_info << e << " " << elemental_error[e] << "\n";
1191  }
1192 
1193  //Now adapt the mesh
1194  Time_mesh_pt->adapt(Problem_pt->communicator_pt(),elemental_error);
1195 
1196  //I seem to have remunged the data,
1197  //Now let's pretend that we have done the the error adaptation
1198  //Time_mesh_pt->refine_uniformly();
1199 
1200  //Let's just refine the central elements twice
1201  //Vector<unsigned> refine_elements;
1202  //refine_elements.push_back(0);
1203  //refine_elements.push_back(9);
1204  //Time_mesh_pt->refine_selected_elements(refine_elements);
1205  //refine_elements.clear();
1206  //refine_elements.push_back(0);
1207  //refine_elements.push_back(1);
1208  //refine_elements.push_back(10);
1209  //refine_elements.push_back(11);
1210  //Time_mesh_pt->refine_selected_elements(refine_elements);
1211 
1212  /* std::ofstream munge("data_refine.dat");
1213  const unsigned n_time_element = Time_mesh_pt->nelement();
1214  for(unsigned e=0;e<n_time_element;e++)
1215  {
1216  const unsigned n_node = Time_mesh_pt->nnode();
1217  for(unsigned n=0;n<n_node;n++)
1218  {
1219  Node* const nod_pt = Time_mesh_pt->node_pt(n);
1220  munge << nod_pt->x(0) << " ";
1221  const unsigned n_space_storage = nod_pt->ntstorage();
1222  for(unsigned t=0;t<n_space_storage;t++)
1223  {
1224  munge << nod_pt->value(t,0) << " ";
1225  }
1226  munge << std::endl;
1227  }
1228  }
1229  munge.close();*/
1230 
1231  //Now we need to put the refined data back into the problem
1232 
1233  //Now we need to number the mesh
1234  //Dummy dof pointer
1235  Vector<double*> dummy_dof_pt;
1236  Time_mesh_pt->assign_global_eqn_numbers(dummy_dof_pt);
1237  //Assign the local equation numbers
1238  Time_mesh_pt->assign_local_eqn_numbers(false);
1239 
1240  //Find the number of temporal degrees of freedom
1241  N_tstorage = dummy_dof_pt.size();
1242  //and new number of elements
1243  N_element_in_period = Time_mesh_pt->nelement();
1244 
1245  //Create the new "fake" timestepper
1246  PeriodicOrbitTimeDiscretisation* periodic_time_stepper_pt =
1248 
1249  //Loop over the temporal elements and set the pointers
1250  for(unsigned e=0;e<N_element_in_period;e++)
1251  {
1254  (Time_mesh_pt->element_pt(e));
1255 
1256  //Set the time and the timestepper
1257  el_pt->time_pt() = Problem_pt->time_pt();
1258  el_pt->time_stepper_pt() = periodic_time_stepper_pt;
1259 
1260  //Set the number of temporal degrees of freedom
1261  el_pt->set_ntstorage(N_tstorage);
1262  //Set the frequency
1263  el_pt->omega_pt() = &Omega;
1264  }
1265 
1266  //We now need to do something much more drastic which is to loop over all
1267  //our the data in the problem and change the timestepper, which is going
1268  //to be a real pain when I start to worry about halo nodes, etc.
1269 
1270  //Will need to use the appropriate mesh-level functions that have
1271  //not been written yet ..
1272 
1273  //Let's just break if there are submeshes
1274  if(Problem_pt->nsub_mesh() > 0)
1275  {
1276  throw OomphLibError(
1277  "PeriodicOrbitHandler can't cope with submeshes yet",
1278  OOMPH_CURRENT_FUNCTION,
1279  OOMPH_EXCEPTION_LOCATION);
1280  }
1281 
1282  //OK now we have only one mesh
1283  for(unsigned n=0;n<n_node;n++)
1284  {
1285  Node* const nod_pt = Problem_pt->mesh_pt()->node_pt(n);
1286  nod_pt->set_time_stepper(periodic_time_stepper_pt,false);
1287  }
1288 
1289  for(unsigned e=0;e<n_space_element;e++)
1290  {
1291  unsigned n_internal =
1293  for(unsigned i=0;i<n_internal;i++)
1294  {
1295  //Cache the internal data
1296  Data* const data_pt = Problem_pt->mesh_pt()->element_pt(e)
1297  ->internal_data_pt(i);
1298  //and set the timestepper
1299  data_pt->set_time_stepper(periodic_time_stepper_pt,false);
1300  }
1301  }
1302 
1303  //Now I can delete the old timestepper and switch
1304  delete Time_stepper_pt;
1305  Time_stepper_pt = periodic_time_stepper_pt;
1306 
1307  //Need to reassign equation numbers so that the DOF pointer, points to the
1308  //newly allocated storage
1309  oomph_info << "Re-allocated " << Problem_pt->assign_eqn_numbers()
1310  << " equation numbers\n";
1311 
1312  //Now's let's add all the unknowns to the problem
1313  Problem_pt->Dof_pt.resize(Ndof*N_tstorage + 1);
1314  //This is reasonably straight forward using pointer arithmetic
1315  //but this does rely on knowing how the data is stored in the
1316  //Nodes which is a little nasty
1317  for(unsigned i=0;i<N_tstorage;i++)
1318  {
1319  unsigned offset = Ndof*i;
1320  for(unsigned n=0;n<Ndof;n++)
1321  {
1322  Problem_pt->Dof_pt[offset + n] = Problem_pt->Dof_pt[n] + i;
1323  }
1324  }
1325 
1326  //Add the frequency of the orbit to the unknowns
1328 
1329  //Rebuild everything
1331  Ndof*N_tstorage+1,true);
1332 
1333 
1334  //Now finally transfer the solution accross
1335 
1336  //Now I "just" copy the values into the new storage
1337  count=0;
1338  for(unsigned i=0;i<n_global_data;i++)
1339  {
1340  Data* const glob_data_pt = Problem_pt->global_data_pt(i);
1341  const unsigned n_value = glob_data_pt->nvalue();
1342  for(unsigned j=0;j<n_value;j++)
1343  {
1344  for(unsigned t=0;t<N_tstorage;t++)
1345  {
1346  //Some heavy assumptions here about the time mesh, but that's OK
1347  //because I know exactly how it's laid out
1348  glob_data_pt->set_value(t,j,
1349  Time_mesh_pt->node_pt(t)->value(count,0));
1350  }
1351  ++count;
1352  }
1353  }
1354 
1355  //Now the nodal data
1356  for(unsigned n=0;n<n_node;n++)
1357  {
1358  Node* const nod_pt = Problem_pt->mesh_pt()->node_pt(n);
1359  const unsigned n_value = nod_pt->nvalue();
1360  for(unsigned i=0;i<n_value;i++)
1361  {
1362  for(unsigned t=0;t<N_tstorage;t++)
1363  {
1364  nod_pt->set_value(t,i,Time_mesh_pt->node_pt(t)->value(count,0));
1365  }
1366  ++count;
1367  }
1368 
1369  //Now deal with the solid node data
1370  SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
1371  if(solid_nod_pt!=0)
1372  {
1373  const unsigned n_solid_value =
1374  solid_nod_pt->variable_position_pt()->nvalue();
1375  for(unsigned i=0;i<n_solid_value;i++)
1376  {
1377  for(unsigned t=0;t<N_tstorage;t++)
1378  {
1379  solid_nod_pt->variable_position_pt()->
1380  set_value(t,i,
1381  Time_mesh_pt->node_pt(t)->value(count,0));
1382  }
1383  ++count;
1384  }
1385  }
1386  }
1387 
1388  //Now just do the internal data
1389  for(unsigned e=0;e<n_space_element;e++)
1390  {
1391  const unsigned n_internal =
1393  for(unsigned i=0;i<n_internal;i++)
1394  {
1395  Data* const internal_dat_pt =
1397 
1398  const unsigned n_value = internal_dat_pt->nvalue();
1399  for(unsigned j=0;j<n_value;j++)
1400  {
1401  for(unsigned t=0;t<N_tstorage;t++)
1402  {
1403  internal_dat_pt->
1404  set_value(t,j,
1405  Time_mesh_pt->node_pt(t)->value(count,0));
1406  }
1407  ++count;
1408  }
1409  }
1410  }
1411 
1412 
1413  //Now I should be able to delete the fake time timestepper
1414  n_time_node = Time_mesh_pt->nnode();
1415  for(unsigned t=0;t<n_time_node;t++)
1416  {
1417  Time_mesh_pt->node_pt(t)->set_time_stepper(
1418  &Mesh::Default_TimeStepper,false);
1419  }
1420  //Delete the fake timestepper
1421  delete fake_space_time_stepper_pt;
1422 
1423  //Set the initial unkowns to be the original problem
1424  Previous_dofs.resize(Ndof*N_tstorage+1,0.0);
1426  }
1427 
1428 
1429  /// \short Destructor, destroy the time mesh
1431 };
1432 
1433 
1435 {
1436 public:
1437 
1439 
1440 
1441 
1442  ///Interface to get the current value of all (internal and shared) unknowns
1444 
1445  ///Interface to get the current value of the time derivative of
1446  /// all (internal and shared) unknowns
1448 
1449 
1450  ///Get the inner product matrix
1451  virtual void get_inner_product_matrix(DenseMatrix<double> &inner_product)
1452  {
1453  const unsigned n_dof = this->ndof();
1454  inner_product.initialise(0.0);
1455  for(unsigned i=0;i<n_dof;i++) {inner_product(i,i) = 1.0;}
1456  }
1457 
1458  virtual void spacetime_output(std::ostream &outilfe, const unsigned &Nplot,
1459  const double &time=0.0)
1460  {}
1461 
1462 
1463 };
1464 
1465 
1466 
1467 
1468 
1469 
1470 }
1471 
1472 #endif
virtual void get_inner_product_matrix(DenseMatrix< double > &inner_product)
Get the inner product matrix.
void assign_initial_positions_impulsive(Node *const &node_pt)
Broken initialisation of the positions for the node corresponding to an impulsive start...
double(* InitialConditionFctPt)(const double &t)
Typedef for function that returns the (scalar) initial value at a given value of the continuous time ...
void adapt_temporal_mesh()
Adapt the time mesh.
A Generalised Element class.
Definition: elements.h:76
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &value)
Return the temporal dummy values.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void fill_in_contribution_to_integrated_residuals(PeriodicOrbitAssemblyHandlerBase *const &assembly_handler_pt, GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
void get_previous_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > &dofs)
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
Problem * Problem_pt
Pointer to the problem.
void initialise(const T &val)
Initialize all values in the matrix to val.
Definition: matrices.h:523
void set_previous_dofs_to_current_dofs()
Update the previous dofs.
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:497
Mesh * Basic_time_mesh_pt
Storage for the mesh of temporal elements with a simple mesh pointer.
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation...
Definition: nodes.cc:2314
PeriodicOrbitTimeDiscretisation(const unsigned &n_tstorage)
Constructor for the case when we allow adaptive timestepping.
bool is_pinned(const unsigned &i) const
Test whether the i-th variable is pinned (1: true; 0: false).
Definition: nodes.h:403
unsigned num_Z2_flux_terms()
Number of flux terms for Z2 error estimation This will be used to represent all spatial values...
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors. ...
Definition: mesh.h:85
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
unsigned N_element_in_period
Storage for number of elements in the period.
Time *const & time_pt() const
Return the pointer to the global time (const version)
double omega()
Return the frequency.
unsigned ndt() const
Number of timestep increments that need to be stored by the scheme.
Refineable version of the OneDMesh.
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:847
PeriodicOrbitTimeDiscretisation * Time_stepper_pt
Pointer to the timestepper.
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
cstr elem_len * i
Definition: cfortran.h:607
double & time()
Return current value of continous time.
Definition: timesteppers.h:324
virtual void get_non_external_dofs(Vector< double > &u)
Interface to get the current value of all (internal and shared) unknowns.
Class to keep track of discrete/continous time. It is essential to have a single Time object when usi...
Definition: timesteppers.h:67
std::string Type
String that indicates the type of the timestepper (e.g. "BDF", "Newmark", etc.)
Definition: timesteppers.h:231
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
const double Pi
50 digits from maple
void output(std::ostream &outfile)
Function to return the number of values.
void operator=(const PeriodicOrbitTimeDiscretisation &)
Broken assignment operator.
The Problem class.
Definition: problem.h:152
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
LinearAlgebraDistribution * Dof_distribution_pt
The distribution of the DOFs in this problem. This object is created in the Problem constructor and s...
Definition: problem.h:457
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
double time() const
Return the global time, accessed via the time pointer.
A general Finite Element class.
Definition: elements.h:1271
void build(const OomphCommunicator *const comm_pt, const unsigned &first_row, const unsigned &nrow_local, const unsigned &nrow=0)
Sets the distribution. Takes first_row, nrow_local and nrow as arguments. If nrow is not provided or ...
char t
Definition: cfortran.h:572
virtual double dshape_and_dtest_eulerian_at_knot_orbit(const unsigned &ipt, Shape &psi, DShape &dpsidt, Shape &test, DShape &dtestdt) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
unsigned order() const
Return the actual order of the scheme.
Vector< double > Previous_dofs
Storage for the previous solution.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void fill_in_contribution_to_integrated_jacobian(PeriodicOrbitAssemblyHandlerBase *const &assembly_handler_pt, GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1264
OomphInfo oomph_info
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
e
Definition: cfortran.h:575
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1221
double Omega
Storage for the frequency of the orbit (scaled by 2pi)
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get the fluxes for the recovert.
double dshape_and_dtest_eulerian_orbit(const Vector< double > &s, Shape &psi, DShape &dpsidt, Shape &test, DShape &dtestdt) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
Time *& time_pt()
Retun the pointer to the global time.
void orbit_output(GeneralisedElement *const &elem_pt, std::ostream &outfile, const unsigned &n_plot)
unsigned nvertex_node() const
Number of vertex nodes in the element.
double *& omega_pt()
Broken assignment operator.
void shift_time_values(Data *const &data_pt)
Broken shifting of time values.
virtual void get_non_external_ddofs_dt(Vector< double > &du_dt)
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
Definition: nodes.h:1654
double dshape_and_dtest_eulerian_at_knot_orbit(const unsigned &ipt, Shape &psi, DShape &dpsidt, Shape &test, DShape &dtestdt) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt. Return Jacobian.
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
unsigned ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
DenseMatrix< double > Weight
Storage for the weights associated with the timestepper.
Definition: timesteppers.h:227
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Return the contribution to the residuals of the element elem_pt.
unsigned nprev_values() const
Number of previous values available.
double & dof(const unsigned &i)
i-th dof in the problem
Definition: problem.h:1696
PeriodicOrbitEquations(const PeriodicOrbitEquations &dummy)
Broken copy constructor.
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
void discrete_times(Vector< double > &t)
Tell me the times at which you want the solution.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values (1)
unsigned N_tstorage
Storage for the number of unknown time values.
PeriodicOrbitTemporalMesh< SpectralPeriodicOrbitElement< NNODE_1D > > * Time_mesh_pt
Storage for mesh of temporal elements.
static char t char * s
Definition: cfortran.h:572
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
Definition: mesh.h:470
void set_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set a new timestepper by resizing the appropriate storage. If already assigned the equation numbering...
Definition: nodes.cc:392
Time * Time_pt
Pointer to discrete time storage scheme.
Definition: timesteppers.h:224
void assemble_residuals_and_jacobian(PeriodicOrbitAssemblyHandlerBase *const &assembly_handler_pt, GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1446
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
virtual void get_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > &dofs)=0
void assign_initial_values_impulsive(Data *const &data_pt)
Broken initialisation the time-history for the Data values corresponding to an impulsive start...
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
A special temporal mesh class.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
PeriodicOrbitTemporalMesh(const unsigned &n_element)
Constructor, create a 1D mesh from 0 to 1 that is periodic.
unsigned nglobal_data() const
Return the number of global data values.
Definition: problem.h:1581
void get_interpolated_values(const Vector< double > &s, Vector< double > &value)
Return the dummy values.
virtual void set_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > const &dofs)=0
A class that is used to define the functions used to assemble the elemental contributions to the resi...
void initialise(const _Tp &__value)
Iterate over all values and set to the desired value.
Definition: Vector.h:171
SpectralPeriodicOrbitElement()
Constructor: Call constructors for QElement and Poisson equations.
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
void set_ntstorage(const unsigned &n_tstorage)
Set the total number of time storage values.
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:456
virtual void spacetime_output(std::ostream &outilfe, const unsigned &Nplot, const double &time=0.0)
Time * Time_pt
Pointer to global time.
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:828
void get_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > &dofs)
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
Definition: nodes.h:1568
void shape(const Vector< double > &s, Shape &psi) const
Calculate the geometric shape functions at local coordinate s.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Return the global equation number of the local unknown ieqn_local in elem_pt.
Timestepper used to calculate periodic orbits directly. It's not really a "timestepper" per se...
unsigned nsub_mesh() const
Return number of submeshes.
Definition: problem.h:1289
Data *& global_data_pt(const unsigned &i)
Return a pointer to the the i-th global data object.
Definition: problem.h:1557
void orbit_output(GeneralisedElement *const &elem_pt, std::ostream &outfile, const unsigned &n_plot)
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void set_timestepper_weights(const Shape &psi, const DShape &dpsidt)
Set the timestepper weights.
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
Definition: problem.cc:1965
Vector< double * > Dof_pt
Vector of pointers to dofs.
Definition: problem.h:551
SpectralPeriodicOrbitElement(const SpectralPeriodicOrbitElement< NNODE_1D > &dummy)
Broken copy constructor.
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
virtual double dshape_and_dtest_eulerian_orbit(const Vector< double > &s, Shape &psi, DShape &dpsidt, Shape &test, DShape &dtestdt) const =0
Shape/test functions and derivs w.r.t. to global coords at local coord. s; return Jacobian of mapping...
void assign_initial_data_values(Data *const &data_pt, Vector< InitialConditionFctPt > initial_value_fct)
Initialise the time-history for the Data values, corresponding to given time history, specified by Vector of function pointers.
unsigned nrecovery_order()
Order of recovery shape functions.
void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable" for elem_pt.
void orbit_output(std::ostream &outfile, const unsigned &n_plot)
Calculate all desired vectors and matrices provided by the element elem_pt.
void fill_in_generic_residual_contribution_orbit(PeriodicOrbitAssemblyHandlerBase *const &assembly_handler_pt, GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
The routine that actually does all the work!
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
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
~PeriodicOrbitAssemblyHandler()
Destructor, destroy the time mesh.
virtual void get_previous_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > &dofs)=0
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
unsigned Ndof
Store number of degrees of freedom in the original problem.
A general mesh class.
Definition: mesh.h:74
void assemble_residuals(PeriodicOrbitAssemblyHandlerBase *const &assembly_handler_pt, GeneralisedElement *const &elem_pt, Vector< double > &residuals)
void set_dofs_for_element(GeneralisedElement *const elem_pt, Vector< double > const &dofs)
void shift_time_positions(Node *const &node_pt)
Broken shifting of time positions.
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
PeriodicOrbitTimeDiscretisation(const PeriodicOrbitTimeDiscretisation &)
Broken copy constructor.