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-2022 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
44namespace 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 {
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.
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
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:
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,
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.
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
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();
775
776 // Set the number of temporal degrees of freedom
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 {
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
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:1313
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:2593
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:3050
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:1432
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
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:3298
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
A Generalised Element class.
Definition: elements.h:73
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
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:704
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:823
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
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
Definition: mesh.h:460
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:493
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)
Time *const & time_pt() const
Return the pointer to the global time (const version)
Time *& time_pt()
Retun the pointer to the global time.
double *& omega_pt()
Broken assignment operator.
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)
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!
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:151
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:1989
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1674
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1246
double & dof(const unsigned &i)
i-th dof in the problem
Definition: problem.h:1813
unsigned nglobal_data() const
Return the number of global data values.
Definition: problem.h:1686
Data *& global_data_pt(const unsigned &i)
Return a pointer to the the i-th global data object.
Definition: problem.h:1647
Vector< double * > Dof_pt
Vector of pointers to dofs.
Definition: problem.h:554
unsigned nsub_mesh() const
Return number of submeshes.
Definition: problem.h:1323
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:460
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1504
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.
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
Number of vertex nodes in the element.
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,...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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...