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