assembly_handler.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //A class that is used to assemble the linear systems solved
31 //by oomph-lib.
32 
33 //Include guards to prevent multiple inclusion of this header
34 #ifndef OOMPH_ASSEMBLY_HANDLER_CLASS_HEADER
35 #define OOMPH_ASSEMBLY_HANDLER_CLASS_HEADER
36 
37 // Config header generated by autoconfig
38 #ifdef HAVE_CONFIG_H
39  #include <oomph-lib-config.h>
40 #endif
41 
42 //OOMPH-LIB headers
43 #include "matrices.h"
44 #include "linear_solver.h"
46 
47 namespace oomph
48 {
49 //Forward class definition of the element
50  class GeneralisedElement;
51 
52  //Forward class definition of the problem
53  class Problem;
54 
55 //=============================================================
56 /// \short A class that is used to define the functions used to
57 /// assemble the elemental contributions to the
58 /// residuals vector and Jacobian matrix that define
59 /// the problem being solved.
60 /// The main use
61 /// of this class is to assemble and solve the augmented systems
62 /// used in bifurcation detection and tracking. The default
63 /// implementation merely calls the underlying elemental
64 /// functions with no augmentation.
65 //===============================================================
67 {
68  public:
69 
70  ///Empty constructor
72 
73  ///Return the number of degrees of freedom in the element elem_pt
74  virtual unsigned ndof(GeneralisedElement* const &elem_pt);
75 
76  ///\short Return vector of dofs at time level t in the element elem_pt
77  virtual void dof_vector(GeneralisedElement* const &elem_pt,
78  const unsigned &t, Vector<double> &dof);
79 
80  ///\short Return vector of pointers to dofs in the element elem_pt
81  virtual void dof_pt_vector(GeneralisedElement* const &elem_pt,
82  Vector<double*> &dof_pt);
83 
84  /// \short Return the t-th level of storage associated with the i-th
85  /// (local) dof stored in the problem
86  virtual double &local_problem_dof(Problem* const &problem_pt,
87  const unsigned &t,
88  const unsigned &i);
89 
90  /// \short Return the global equation number of the local unknown ieqn_local
91  ///in elem_pt.
92  virtual unsigned long eqn_number(GeneralisedElement* const &elem_pt,
93  const unsigned &ieqn_local);
94 
95  ///Return the contribution to the residuals of the element elem_pt
96  virtual void get_residuals(GeneralisedElement* const &elem_pt,
97  Vector<double> &residuals);
98 
99  /// \short Calculate the elemental Jacobian matrix "d equation
100  /// / d variable" for elem_pt.
101  virtual void get_jacobian(GeneralisedElement* const &elem_pt,
102  Vector<double> &residuals,
103  DenseMatrix<double> &jacobian);
104 
105  /// \short Calculate all desired vectors and matrices
106  /// provided by the element elem_pt.
107  virtual void get_all_vectors_and_matrices(
108  GeneralisedElement* const &elem_pt,
110 
111  /// \short Calculate the derivative of the residuals with respect to
112  /// a parameter
113  virtual void get_dresiduals_dparameter(GeneralisedElement* const &elem_pt,
114  double* const &parameter_pt,
115  Vector<double> &dres_dparam);
116 
117  /// \short Calculate the derivative of the residuals and jacobian
118  /// with respect to a parameter
119  virtual void get_djacobian_dparameter(GeneralisedElement* const &elem_pt,
120  double* const &parameter_pt,
121  Vector<double> &dres_dparam,
122  DenseMatrix<double> &djac_dparam);
123 
124  /// \short Calculate the product of the Hessian (derivative of Jacobian with
125  /// respect to all variables) an eigenvector, Y, and
126  /// other specified vectors, C
127  /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}
128  virtual void get_hessian_vector_products(GeneralisedElement* const &elem_pt,
129  Vector<double> const &Y,
130  DenseMatrix<double> const &C,
131  DenseMatrix<double> &product);
132 
133 
134  /// \short Return an unsigned integer to indicate whether the
135  /// handler is a bifurcation tracking handler. The default
136  /// is zero (not)
137  virtual int bifurcation_type() const {return 0;}
138 
139  /// \short Return a pointer to the
140  /// bifurcation parameter in bifurcation tracking problems
141  virtual double* bifurcation_parameter_pt() const;
142 
143  /// \short Return the eigenfunction(s) associated with the bifurcation that
144  /// has been detected in bifurcation tracking problems
145  virtual void get_eigenfunction(Vector<DoubleVector> &eigenfunction);
146 
147  /// \short Compute the inner products of the given vector of pairs of
148  /// history values over the element.
149  virtual void get_inner_products(GeneralisedElement* const &elem_pt,
150  Vector<std::pair<unsigned,unsigned> >
151  const &history_index,
152  Vector<double> &inner_product);
153 
154  /// \short Compute the vectors that when taken as a dot product with
155  /// other history values give the inner product over the element
156  virtual void get_inner_product_vectors(GeneralisedElement* const &elem_pt,
157  Vector<unsigned> const &history_index,
159  &inner_product_vector);
160 
161 #ifdef OOMPH_HAS_MPI
162 
163  /// \short Function that is used to perform any synchronisation
164  /// required during the solution
165  virtual void synchronise() {}
166 #endif
167 
168  /// \short Empty virtual destructor
169  virtual ~AssemblyHandler() {}
170 };
171 
172 
173 
174 //=============================================================
175 /// \short A class that is used to define the functions used to
176 /// assemble and invert the mass matrix when taking an explicit
177 /// timestep. The idea is simply to replace the jacobian matrix
178 /// with the mass matrix and then our standard linear solvers
179 /// will solve the required system
180 //===============================================================
182 {
183  public:
184 
185  ///Empty Constructor
187 
188  ///Return the number of degrees of freedom in the element elem_pt
189  unsigned ndof(GeneralisedElement* const &elem_pt);
190 
191  ///\short Return the global equation number of the local unknown ieqn_local
192  ///in elem_pt.
193  unsigned long eqn_number(GeneralisedElement* const &elem_pt,
194  const unsigned &ieqn_local);
195 
196  ///\short Return the contribution to the residuals of the element elem_pt
197  ///This is deliberately broken in our eigenproblem
198  void get_residuals(GeneralisedElement* const &elem_pt,
199  Vector<double> &residuals);
200 
201  /// \short Calculate the elemental Jacobian matrix "d equation
202  /// / d variable" for elem_pt. Again deliberately broken in the eigenproblem
203  void get_jacobian(GeneralisedElement* const &elem_pt,
204  Vector<double> &residuals,
205  DenseMatrix<double> &jacobian);
206 
207  /// \short Calculate all desired vectors and matrices
208  /// provided by the element elem_pt.
210  GeneralisedElement* const &elem_pt,
212 
213  /// \short Empty virtual destructor
215 
216 };
217 
218 
219 
220 //=============================================================
221 /// \short A class that is used to define the functions used to
222 /// assemble the elemental contributions to the
223 /// mass matrix and jacobian (stiffness) matrix that define
224 /// a generalised eigenproblem.
225 //===============================================================
227 {
228  /// Storage for the real shift
229  double Sigma_real;
230 
231  public:
232 
233  ///Constructor, sets the value of the real shift
234  EigenProblemHandler(const double &sigma_real) : Sigma_real(sigma_real) {}
235 
236  ///Return the number of degrees of freedom in the element elem_pt
237  unsigned ndof(GeneralisedElement* const &elem_pt);
238 
239  ///\short Return the global equation number of the local unknown ieqn_local
240  ///in elem_pt.
241  unsigned long eqn_number(GeneralisedElement* const &elem_pt,
242  const unsigned &ieqn_local);
243 
244  ///\short Return the contribution to the residuals of the element elem_pt
245  ///This is deliberately broken in our eigenproblem
246  void get_residuals(GeneralisedElement* const &elem_pt,
247  Vector<double> &residuals);
248 
249  /// \short Calculate the elemental Jacobian matrix "d equation
250  /// / d variable" for elem_pt. Again deliberately broken in the eigenproblem
251  void get_jacobian(GeneralisedElement* const &elem_pt,
252  Vector<double> &residuals,
253  DenseMatrix<double> &jacobian);
254 
255  /// \short Calculate all desired vectors and matrices
256  /// provided by the element elem_pt.
258  GeneralisedElement* const &elem_pt,
260 
261  /// \short Empty virtual destructor
263 
264 };
265 
266 
267 //=============================================================
268 /// \short A class that is used to assemble the residuals in
269 /// parallel by overloading the get_all_vectors_and_matrices,
270 /// so that only the residuals are returned. This ensures that
271 /// the (moderately complex) distributed parallel assembly
272 /// loops are only in one place.
273 //===============================================================
275 {
276  ///The original assembly handler
278 
279  public:
280 
281  ///Constructor, set the original assembly handler
282  ParallelResidualsHandler(AssemblyHandler* const &assembly_handler_pt) :
283  Assembly_handler_pt(assembly_handler_pt) {}
284 
285  ///\short Use underlying assembly handler to return the number of
286  /// degrees of freedom in the element elem_pt
287  unsigned ndof(GeneralisedElement* const &elem_pt)
288  {return Assembly_handler_pt->ndof(elem_pt);}
289 
290  /// \short Use underlying AssemblyHandler to return the
291  /// global equation number of the local unknown ieqn_local in elem_pt.
292  unsigned long eqn_number(GeneralisedElement* const &elem_pt,
293  const unsigned &ieqn_local)
294  {return Assembly_handler_pt->eqn_number(elem_pt,ieqn_local);}
295 
296  ///\short Use underlying AssemblyHandler to return
297  /// the contribution to the residuals of the element elem_pt
298  void get_residuals(GeneralisedElement* const &elem_pt,
299  Vector<double> &residuals)
300  {Assembly_handler_pt->get_residuals(elem_pt,residuals);}
301 
302 
303  /// \short Use underlying AssemblyHandler to
304  /// Calculate the elemental Jacobian matrix "d equation
305  /// / d variable" for elem_pt.
306  void get_jacobian(GeneralisedElement* const &elem_pt,
307  Vector<double> &residuals,
308  DenseMatrix<double> &jacobian)
309  {Assembly_handler_pt->get_jacobian(elem_pt,residuals,jacobian);}
310 
311 
312  /// \short Calculate all desired vectors and matrices
313  /// provided by the element elem_pt
314  /// This function calls only the get_residuals function associated
315  /// with the original assembly handler
317  GeneralisedElement* const &elem_pt,
319  {Assembly_handler_pt->get_residuals(elem_pt,vec[0]);}
320 
321  /// \short Empty virtual destructor
323 
324 };
325 
326 
327 
328 //==========================================================================
329 /// \short A class that is used to define the functions used when assembling
330 /// the derivatives of the residuals with respect to a parameter.
331 /// The idea is to replace get_residuals with get_dresiduals_dparameter with
332 /// a particular parameter and assembly handler that are passed on
333 /// assembly.
334 //==========================================================================
336 {
337  ///The value of the parameter
338  double *Parameter_pt;
339 
340  ///The original assembly handler
342 
343  public:
344 
345  ///Store the original assembly handler and parameter
346  ParameterDerivativeHandler(AssemblyHandler* const &assembly_handler_pt,
347  double* const &parameter_pt) :
348  Parameter_pt(parameter_pt), Assembly_handler_pt(assembly_handler_pt)
349  { }
350 
351  ///Return the number of degrees of freedom in the element elem_pt
352  ///Pass through to the original assembly handler
353  unsigned ndof(GeneralisedElement* const &elem_pt)
354  {return Assembly_handler_pt->ndof(elem_pt);}
355 
356  ///\short Return the global equation number of the local unknown ieqn_local
357  ///in elem_pt.Pass through to the original assembly handler
358  unsigned long eqn_number(GeneralisedElement* const &elem_pt,
359  const unsigned &ieqn_local)
360  {return Assembly_handler_pt->eqn_number(elem_pt,ieqn_local);}
361 
362  ///\short Return the contribution to the residuals of the element elem_pt
363  /// by using the derivatives
364  void get_residuals(GeneralisedElement* const &elem_pt,
365  Vector<double> &residuals)
367  residuals);}
368 
369 
370  /// \short Calculate the elemental Jacobian matrix "d equation
371  /// / d variable" for elem_pt.
372  /// Overloaded to return the derivatives wrt the parameter
373  void get_jacobian(GeneralisedElement* const &elem_pt,
374  Vector<double> &residuals,
375  DenseMatrix<double> &jacobian)
377  residuals,jacobian);}
378 
379 };
380 
381 
382 
383 
384 //========================================================================
385 /// A custom linear solver class that is used to solve a block-factorised
386 /// version of the Fold bifurcation detection problem.
387 //========================================================================
389 {
390  ///Pointer to the original linear solver
392 
393  //Pointer to the problem, used in the resolve
395 
396  ///Pointer to the storage for the vector alpha
398 
399  ///Pointer to the storage for the vector e
401 
402 public:
403 
404  ///Constructor, inherits the original linear solver
406  : Linear_solver_pt(linear_solver_pt), Problem_pt(0), Alpha_pt(0), E_pt(0) {}
407 
408  ///Destructor: clean up the allocated memory
410 
411  /// The solve function uses the block factorisation
412  void solve(Problem* const &problem_pt, DoubleVector &result);
413 
414  /// \short The linear-algebra-type solver does not make sense.
415  /// The interface is deliberately broken
416  void solve(DoubleMatrixBase* const &matrix_pt,
417  const DoubleVector &rhs,
418  DoubleVector &result)
419  {
420  throw OomphLibError(
421  "Linear-algebra interface does not make sense for this linear solver\n",
422  OOMPH_CURRENT_FUNCTION,
423  OOMPH_EXCEPTION_LOCATION);
424  }
425 
426  /// \short The linear-algebra-type solver does not make sense.
427  /// The interface is deliberately broken
428  void solve(DoubleMatrixBase* const &matrix_pt,
429  const Vector<double> &rhs,
430  Vector<double> &result)
431  {
432  throw OomphLibError(
433  "Linear-algebra interface does not make sense for this linear solver\n",
434  OOMPH_CURRENT_FUNCTION,
435  OOMPH_EXCEPTION_LOCATION);
436  }
437 
438  /// The resolve function also uses the block factorisation
439  void resolve(const DoubleVector &rhs, DoubleVector &result);
440 
441  /// Access function to the original linear solver
443 
444 };
445 
446 
447 //===================================================================
448 /// A class that is used to assemble the augmented system that defines
449 /// a fold (saddle-node) or limit point. The "standard" problem must
450 /// be a function of a global paramter \f$\lambda\f$, and a
451 /// solution is \f$R(u,\lambda) = 0 \f$ , where \f$ u \f$ are the unknowns
452 /// in the problem. A limit point is formally specified by the augmented
453 /// system of size \f$ 2N+1 \f$
454 /// \f[ R(u,\lambda) = 0, \f]
455 /// \f[ Jy = 0, \f]
456 /// \f[\phi\cdot y = 1. \f]
457 /// In the above \f$ J \f$ is the usual Jacobian matrix, \f$ dR/du \f$, and
458 /// \f$ \phi \f$ is a constant vector that is chosen to
459 /// ensure that the null vector, \f$ y \f$, is not trivial.
460 //======================================================================
462  {
464 
465  /// A little private enum to determine whether we are solving
466  /// the block system or not
468 
469  /// \short Integer flag to indicate which system should be assembled.
470  /// There are three possibilities. The full augmented system (0),
471  /// the non-augmented jacobian system (1), a system in which the
472  /// jacobian is augmented by 1 row and column to ensure that it is
473  /// non-singular (2). See the enum above
475 
476  ///Pointer to the problem
478 
479  /// \short Store the number of degrees of freedom in the non-augmented
480  ///problem
481  unsigned Ndof;
482 
483  /// \short A constant vector used to ensure that the null vector
484  /// is not trivial
486 
487  /// \short Storage for the null vector
489 
490  /// \short A vector that is used to determine how many elements
491  /// contribute to a particular equation. It is used to ensure
492  /// that the global system is correctly formulated.
494 
495  /// \short Storage for the pointer to the parameter
496  double *Parameter_pt;
497 
498  public:
499 
500  ///\short Constructor:
501  /// initialise the fold handler, by setting initial guesses
502  /// for Y, Phi and calculating count. If the system changes, a new
503  /// fold handler must be constructed
504  FoldHandler(Problem* const &problem_pt,
505  double* const &parameter_pt);
506 
507 
508  /// Constructor in which initial eigenvector can be passed
509  FoldHandler(Problem* const &problem_pt,
510  double* const &parameter_pt,
511  const DoubleVector &eigenvector);
512 
513  /// Constructor in which initial eigenvector
514  /// and normalisation can be passed
515  FoldHandler(Problem* const &problem_pt,
516  double* const &parameter_pt,
517  const DoubleVector &eigenvector,
518  const DoubleVector &normalisation);
519 
520 
521  /// \short Destructor, return the problem to its original state
522  /// before the augmented system was added
523  ~FoldHandler();
524 
525  ///Get the number of elemental degrees of freedom
526  unsigned ndof(GeneralisedElement* const &elem_pt);
527 
528  ///Get the global equation number of the local unknown
529  unsigned long eqn_number(GeneralisedElement* const &elem_pt,
530  const unsigned &ieqn_local);
531 
532  ///Get the residuals
533  void get_residuals(GeneralisedElement* const &elem_pt,
534  Vector<double> &residuals);
535 
536  /// \short Calculate the elemental Jacobian matrix "d equation
537  /// / d variable".
538  void get_jacobian(GeneralisedElement* const &elem_pt,
539  Vector<double> &residuals,
540  DenseMatrix<double> &jacobian);
541 
542  /// \short Overload the derivatives of the residuals with respect to
543  /// a parameter to apply to the augmented system
544  void get_dresiduals_dparameter(GeneralisedElement* const &elem_pt,
545  double* const &parameter_pt,
546  Vector<double> &dres_dparam);
547 
548  /// \short Overload the derivative of the residuals and jacobian
549  /// with respect to a parameter so that it breaks
550  void get_djacobian_dparameter(GeneralisedElement* const &elem_pt,
551  double* const &parameter_pt,
552  Vector<double> &dres_dparam,
553  DenseMatrix<double> &djac_dparam);
554 
555  /// \short Overload the hessian vector product function so that
556  /// it breaks
557  void get_hessian_vector_products(GeneralisedElement* const &elem_pt,
558  Vector<double> const &Y,
559  DenseMatrix<double> const &C,
560  DenseMatrix<double> &product);
561 
562  ///\short Indicate that we are tracking a fold bifurcation by returning 1
563  int bifurcation_type() const {return 1;}
564 
565  /// \short Return a pointer to the
566  /// bifurcation parameter in bifurcation tracking problems
567  double* bifurcation_parameter_pt() const
568  {return Parameter_pt;}
569 
570  /// \short Return the eigenfunction(s) associated with the bifurcation that
571  /// has been detected in bifurcation tracking problems
572  void get_eigenfunction(Vector<DoubleVector> &eigenfunction);
573 
574  /// \short Set to solve the augmented block system
576 
577  /// \short Set to solve the block system
578  void solve_block_system();
579 
580  /// \short Solve non-block system
581  void solve_full_system();
582 
583  };
584 
585 
586 
587 //========================================================================
588 /// A custom linear solver class that is used to solve a block-factorised
589 /// version of the PitchFork bifurcation detection problem.
590 //========================================================================
592 {
593  ///Pointer to the original linear solver
595 
596  ///Pointer to the problem, used in the resolve
598 
599  ///Pointer to the storage for the vector b
601 
602  ///Pointer to the storage for the vector c
604 
605  ///Pointer to the storage for the vector d
607 
608  ///Pointer to the storage for the vector of derivatives with respect
609  ///to the bifurcation parameter
611 
612 public:
613 
614  ///Constructor, inherits the original linear solver
616  : Linear_solver_pt(linear_solver_pt), Problem_pt(0),
617  B_pt(0), C_pt(0), D_pt(0), dJy_dparam_pt(0) {}
618 
619  ///Destructor: clean up the allocated memory
621 
622  /// The solve function uses the block factorisation
623  void solve(Problem* const &problem_pt, DoubleVector &result);
624 
625  /// \short The linear-algebra-type solver does not make sense.
626  /// The interface is deliberately broken
627  void solve(DoubleMatrixBase* const &matrix_pt,
628  const DoubleVector &rhs,
629  DoubleVector &result)
630  {
631  throw OomphLibError(
632  "Linear-algebra interface does not make sense for this linear solver\n",
633  OOMPH_CURRENT_FUNCTION,
634  OOMPH_EXCEPTION_LOCATION);
635  }
636 
637  /// \short The linear-algebra-type solver does not make sense.
638  /// The interface is deliberately broken
639  void solve(DoubleMatrixBase* const &matrix_pt,
640  const Vector<double> &rhs,
641  Vector<double> &result)
642  {
643  throw OomphLibError(
644  "Linear-algebra interface does not make sense for this linear solver\n",
645  OOMPH_CURRENT_FUNCTION,
646  OOMPH_EXCEPTION_LOCATION);
647  }
648 
649 
650  /// The resolve function also uses the block factorisation
651  void resolve(const DoubleVector &rhs, DoubleVector &result);
652 
653  /// Access function to the original linear solver
655 
656 };
657 
658 
659 
660 
661 //========================================================================
662 /// A custom linear solver class that is used to solve a block-factorised
663 /// version of the PitchFork bifurcation detection problem.
664 //========================================================================
666 {
667  ///Pointer to the original linear solver
669 
670  ///Pointer to the problem, used in the resolve
672 
673  ///Pointer to the storage for the vector alpha
675 
676  ///Pointer to the storage for the vector e
678 
679 public:
680 
681  ///Constructor, inherits the original linear solver
683  : Linear_solver_pt(linear_solver_pt), Problem_pt(0),
684  Alpha_pt(0), E_pt(0) {}
685 
686  ///Destructor: clean up the allocated memory
688 
689  /// The solve function uses the block factorisation
690  void solve(Problem* const &problem_pt, DoubleVector &result);
691 
692  /// \short The linear-algebra-type solver does not make sense.
693  /// The interface is deliberately broken
694  void solve(DoubleMatrixBase* const &matrix_pt,
695  const DoubleVector &rhs,
696  DoubleVector &result)
697  {
698  throw OomphLibError(
699  "Linear-algebra interface does not make sense for this linear solver\n",
700  OOMPH_CURRENT_FUNCTION,
701  OOMPH_EXCEPTION_LOCATION);
702  }
703 
704  /// \short The linear-algebra-type solver does not make sense.
705  /// The interface is deliberately broken
706  void solve(DoubleMatrixBase* const &matrix_pt,
707  const Vector<double> &rhs,
708  Vector<double> &result)
709  {
710  throw OomphLibError(
711  "Linear-algebra interface does not make sense for this linear solver\n",
712  OOMPH_CURRENT_FUNCTION,
713  OOMPH_EXCEPTION_LOCATION);
714  }
715 
716 
717  /// The resolve function also uses the block factorisation
718  void resolve(const DoubleVector &rhs, DoubleVector &result);
719 
720  /// Access function to the original linear solver
722 
723 };
724 
725 //========================================================================
726 /// A class that is used to assemble the augmented system that defines
727 /// a pitchfork (symmetry-breaking) bifurcation. The "standard" problem
728 /// must be a function of a global parameter \f$ \lambda \f$ and a solution
729 /// is \f$R(u,\lambda) = 0\f$, where \f$u\f$ are the unknowns in the
730 /// problem. A pitchfork bifurcation may be specified by the augmented
731 /// system of size \f$2N+2\f$.
732 /// \f[ R(u,\lambda) + \sigma \psi = 0,\f]
733 /// \f[ Jy = 0,\f]
734 /// \f[ \langle u, \phi \rangle = 0, \f]
735 /// \f[ \phi \cdot y = 1.\f]
736 /// In the abovem \f$J\f$ is the usual Jacobian matrix, \f$dR/du\f$
737 /// and \f$\phi\f$ is a constant vector that is chosen to ensure that
738 /// the null vector, \f$y\f$, is not trivial.
739 /// Here \f$\sigma \f$ is a slack variable that is used to enforce the
740 /// constraint \f$ \langle u, \phi \rangle = 0 \f$ --- the inner product
741 /// of the solution with the chosen symmetry vector is zero. At the
742 /// bifurcation \f$ \sigma \f$ should be very close to zero. Note that
743 /// this formulation means that any odd symmetry in the problem must
744 /// be about zero. Moreover, we use the dot product of two vectors to
745 /// calculate the inner product, rather than an integral over the
746 /// domain and so the mesh must be symmetric in the appropriate directions.
747 //========================================================================
749  {
752 
753  /// A little private enum to determine whether we are solving
754  /// the block system or not
756 
757  /// \short Integer flag to indicate which system should be assembled.
758  /// There are three possibilities. The full augmented system (0),
759  /// the non-augmented jacobian system (1), a system in which the
760  /// jacobian is augmented by 1 row and column to ensure that it is
761  /// non-singular (2). See the enum above
763 
764  ///Pointer to the problem
766 
767  ///Pointer to the underlying (original) assembly handler
769 
770  /// \short Store the number of degrees of freedom in the non-augmented
771  ///problem
772  unsigned Ndof;
773 
774  /// \short Store the original dof distribution
776 
777  /// \short The augmented distribution
779 
780  /// \short A slack variable used to specify the amount of antisymmetry
781  /// in the solution.
782  double Sigma;
783 
784  /// \short Storage for the null vector
786 
787  /// \short A constant vector that is specifies the symmetry being broken
789 
790  /// \short A constant vector used to ensure that the null vector
791  /// is not trivial
793 
794  /// \short A vector that is used to determine how many elements
795  /// contribute to a particular equation. It is used to ensure
796  /// that the global system is correctly formulated.
797  /// This should really be an integer, but its double so that
798  /// the distribution can be used
800 
801  /// \short A vector that is used to map the global equations to their
802  /// actual location in a distributed problem
804 
805  /// \short Storage for the pointer to the parameter
806  double *Parameter_pt;
807 
808  // \short The total number of elements in the problem
809  unsigned Nelement;
810 
811 #ifdef OOMPH_HAS_MPI
812  /// \short Boolean to indicate whether the problem is distributed
814 #endif
815 
816  /// \short Function that is used to return map the global equations
817  /// using the simplistic numbering scheme into the actual distributed
818  /// scheme
819  inline unsigned global_eqn_number(const unsigned &i)
820  {
821 #ifdef OOMPH_HAS_MPI
822  //If the problem is distributed I have to do something
823  if(Distributed)
824  {
825  return Global_eqn_number[i];
826  }
827  //Otherwise it's just i
828  else
829 #endif
830  {
831  return i;
832  }
833  }
834 
835  public:
836 
837  ///Constructor, initialise the systems
838  PitchForkHandler(Problem* const &problem_pt,
839  AssemblyHandler* const &assembly_handler_pt,
840  double* const &parameter_pt,
841  const DoubleVector &symmetry_vector);
842 
843  /// Destructor, return the problem to its original state,
844  /// before the augmented system was added
846 
847  //Has this been called
848 
849  ///Get the number of elemental degrees of freedom
850  unsigned ndof(GeneralisedElement* const &elem_pt);
851 
852  ///Get the global equation number of the local unknown
853  unsigned long eqn_number(GeneralisedElement* const &elem_pt,
854  const unsigned &ieqn_local);
855 
856  ///Get the residuals
857  void get_residuals(GeneralisedElement* const &elem_pt,
858  Vector<double> &residuals);
859 
860  /// \short Calculate the elemental Jacobian matrix "d equation
861  /// / d variable".
862  void get_jacobian(GeneralisedElement* const &elem_pt,
863  Vector<double> &residuals,
864  DenseMatrix<double> &jacobian);
865 
866  /// \short Overload the derivatives of the residuals with respect to
867  /// a parameter to apply to the augmented system
868  void get_dresiduals_dparameter(GeneralisedElement* const &elem_pt,
869  double* const &parameter_pt,
870  Vector<double> &dres_dparam);
871 
872  /// \short Overload the derivative of the residuals and jacobian
873  /// with respect to a parameter so that it breaks
874  void get_djacobian_dparameter(GeneralisedElement* const &elem_pt,
875  double* const &parameter_pt,
876  Vector<double> &dres_dparam,
877  DenseMatrix<double> &djac_dparam);
878 
879  /// \short Overload the hessian vector product function so that
880  /// it breaks
881  void get_hessian_vector_products(GeneralisedElement* const &elem_pt,
882  Vector<double> const &Y,
883  DenseMatrix<double> const &C,
884  DenseMatrix<double> &product);
885 
886 
887  ///\short Indicate that we are tracking a pitchfork
888  /// bifurcation by returning 2
889  int bifurcation_type() const {return 2;}
890 
891  /// \short Return a pointer to the
892  /// bifurcation parameter in bifurcation tracking problems
893  double* bifurcation_parameter_pt() const
894  {return Parameter_pt;}
895 
896  /// \short Return the eigenfunction(s) associated with the bifurcation that
897  /// has been detected in bifurcation tracking problems
898  void get_eigenfunction(Vector<DoubleVector> &eigenfunction);
899 
900 #ifdef OOMPH_HAS_MPI
901  /// \short Function that is used to perform any synchronisation
902  /// required during the solution
903  void synchronise();
904 #endif
905 
906 
907  /// \short Set to solve the augmented block system
909 
910  /// \short Set to solve the block system
911  void solve_block_system();
912 
913  /// \short Solve non-block system
914  void solve_full_system();
915 
916  };
917 
918 //========================================================================
919 /// A custom linear solver class that is used to solve a block-factorised
920 /// version of the Hopf bifurcation detection problem.
921 //========================================================================
923 {
924  ///Pointer to the original linear solver
926 
927  ///Pointer to the problem, used in the resolve
929 
930  ///Pointer to the storage for the vector a
932 
933  ///Pointer to the storage for the vector e (0 to n-1)
935 
936  ///Pointer to the storage for the vector g (0 to n-1)
938 
939 public:
940 
941  ///Constructor, inherits the original linear solver
943  : Linear_solver_pt(linear_solver_pt), Problem_pt(0),
944  A_pt(0), E_pt(0), G_pt(0) {}
945 
946  ///Destructor: clean up the allocated memory
948 
949  /// Solve for two right hand sides
950  void solve_for_two_rhs(Problem* const &problem_pt,
951  DoubleVector &result,
952  const DoubleVector &rhs2,
953  DoubleVector &result2);
954 
955  /// The solve function uses the block factorisation
956  void solve(Problem* const &problem_pt, DoubleVector &result);
957 
958  /// \short The linear-algebra-type solver does not make sense.
959  /// The interface is deliberately broken
960  void solve(DoubleMatrixBase* const &matrix_pt,
961  const DoubleVector &rhs,
962  DoubleVector &result)
963  {
964  throw OomphLibError(
965  "Linear-algebra interface does not make sense for this linear solver\n",
966  OOMPH_CURRENT_FUNCTION,
967  OOMPH_EXCEPTION_LOCATION);
968  }
969 
970  /// \short The linear-algebra-type solver does not make sense.
971  /// The interface is deliberately broken
972  void solve(DoubleMatrixBase* const &matrix_pt,
973  const Vector<double> &rhs,
974  Vector<double> &result)
975  {
976  throw OomphLibError(
977  "Linear-algebra interface does not make sense for this linear solver\n",
978  OOMPH_CURRENT_FUNCTION,
979  OOMPH_EXCEPTION_LOCATION);
980  }
981 
982 
983  /// The resolve function also uses the block factorisation
984  void resolve(const DoubleVector &rhs, DoubleVector &result);
985 
986  /// Access function to the original linear solver
988 
989 };
990 
991 
992 //===============================================================
993 /// A class that is used to assemble the augmented system that defines
994 /// a Hopf bifurcation. The "standard" problem
995 /// must be a function of a global parameter \f$ \lambda \f$ and a solution
996 /// is \f$ R(u,\lambda) = 0 \f$, where \f$ u \f$ are the unknowns in the
997 /// problem. A Hopf bifurcation may be specified by the augmented
998 /// system of size \f$ 3N+2 \f$.
999 /// \f[ R(u,\lambda) = 0, \f]
1000 /// \f[ J\phi + \omega M \psi = 0, \f]
1001 /// \f[ J\psi - \omega M \phi = 0, \f]
1002 /// \f[ c \cdot \phi = 1. \f]
1003 /// \f[ c \cdot \psi = 0. \f]
1004 /// In the above \f$ J \f$ is the usual Jacobian matrix, \f$ dR/du \f$
1005 /// and \f$ M \f$ is the mass matrix that multiplies the time derivative terms.
1006 /// \f$ \phi + i\psi \f$ is the (complex) null vector of the complex matrix
1007 /// \f$ J - i\omega M \f$, where \f$ \omega \f$ is the critical frequency.
1008 /// \f$ c \f$ is a constant vector that is used to ensure that the null vector
1009 /// is non-trivial.
1010 //===========================================================================
1012  {
1014 
1015  ///\short Integer flag to indicate which system should be assembled.
1016  /// There are three possibilities. The full augmented system (0),
1017  /// the non-augmented jacobian system (1), and complex
1018  /// system (2), where the matrix is a combination of the jacobian
1019  /// and mass matrices.
1021 
1022  ///Pointer to the problem
1024 
1025  ///Pointer to the parameter
1026  double *Parameter_pt;
1027 
1028  /// \short Store the number of degrees of freedom in the non-augmented
1029  ///problem
1030  unsigned Ndof;
1031 
1032  /// \short The critical frequency of the bifurcation
1033  double Omega;
1034 
1035  /// \short The real part of the null vector
1037 
1038  /// \short The imaginary part of the null vector
1040 
1041  /// \short A constant vector used to ensure that the null vector is
1042  /// not trivial
1044 
1045  /// \short A vector that is used to determine how many elements
1046  /// contribute to a particular equation. It is used to ensure
1047  /// that the global system is correctly formulated.
1049 
1050  public:
1051 
1052  /// Constructor
1053  HopfHandler(Problem* const &problem_pt, double* const &parameter_pt);
1054 
1055  /// Constructor with initial guesses for the frequency and null
1056  /// vectors, such as might be provided by an eigensolver
1057  HopfHandler(Problem* const &problem_pt, double* const &paramter_pt,
1058  const double &omega, const DoubleVector &phi,
1059  const DoubleVector &psi);
1060 
1061  /// Destructor, return the problem to its original state,
1062  /// before the augmented system was added
1063  ~HopfHandler();
1064 
1065  ///Get the number of elemental degrees of freedom
1066  unsigned ndof(GeneralisedElement* const &elem_pt);
1067 
1068  ///Get the global equation number of the local unknown
1069  unsigned long eqn_number(GeneralisedElement* const &elem_pt,
1070  const unsigned &ieqn_local);
1071 
1072  ///Get the residuals
1073  void get_residuals(GeneralisedElement* const &elem_pt,
1074  Vector<double> &residuals);
1075 
1076  /// \short Calculate the elemental Jacobian matrix "d equation
1077  /// / d variable".
1078  void get_jacobian(GeneralisedElement* const &elem_pt,
1079  Vector<double> &residuals,
1080  DenseMatrix<double> &jacobian);
1081 
1082  /// \short Overload the derivatives of the residuals with respect to
1083  /// a parameter to apply to the augmented system
1084  void get_dresiduals_dparameter(GeneralisedElement* const &elem_pt,
1085  double* const &parameter_pt,
1086  Vector<double> &dres_dparam);
1087 
1088  /// \short Overload the derivative of the residuals and jacobian
1089  /// with respect to a parameter so that it breaks
1090  void get_djacobian_dparameter(GeneralisedElement* const &elem_pt,
1091  double* const &parameter_pt,
1092  Vector<double> &dres_dparam,
1093  DenseMatrix<double> &djac_dparam);
1094 
1095  /// \short Overload the hessian vector product function so that
1096  /// it breaks
1097  void get_hessian_vector_products(GeneralisedElement* const &elem_pt,
1098  Vector<double> const &Y,
1099  DenseMatrix<double> const &C,
1100  DenseMatrix<double> &product);
1101 
1102  ///\short Indicate that we are tracking a Hopf
1103  /// bifurcation by returning 3
1104  int bifurcation_type() const {return 3;}
1105 
1106  /// \short Return a pointer to the
1107  /// bifurcation parameter in bifurcation tracking problems
1108  double* bifurcation_parameter_pt() const
1109  {return Parameter_pt;}
1110 
1111  /// \short Return the eigenfunction(s) associated with the bifurcation that
1112  /// has been detected in bifurcation tracking problems
1113  void get_eigenfunction(Vector<DoubleVector> &eigenfunction);
1114 
1115  /// \short Return the frequency of the bifurcation
1116  const double &omega() const {return Omega;}
1117 
1118  /// \short Set to solve the standard system
1119  void solve_standard_system();
1120 
1121  /// \short Set to solve the complex system
1122  void solve_complex_system();
1123 
1124  /// \short Solve non-block system
1125  void solve_full_system();
1126 
1127  };
1128 
1129 
1130 
1131 
1132 }
1133 #endif
A Generalised Element class.
Definition: elements.h:76
DoubleVectorWithHaloEntries Count
A vector that is used to determine how many elements contribute to a particular equation. It is used to ensure that the global system is correctly formulated. This should really be an integer, but its double so that the distribution can be used.
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
void get_all_vectors_and_matrices(GeneralisedElement *const &elem_pt, Vector< Vector< double > > &vec, Vector< DenseMatrix< double > > &matrix)
Calculate all desired vectors and matrices provided by the element elem_pt.
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. Overloaded to return ...
void solve_block_system()
Set to solve the block system.
~BlockPitchForkLinearSolver()
Destructor: clean up the allocated memory.
double Sigma
A slack variable used to specify the amount of antisymmetry in the solution.
void get_dresiduals_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam)
Overload the derivatives of the residuals with respect to.
AssemblyHandler * Assembly_handler_pt
Pointer to the underlying (original) assembly handler.
DoubleVector * E_pt
Pointer to the storage for the vector e (0 to n-1)
DoubleVector * A_pt
Pointer to the storage for the vector a.
void get_dresiduals_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam)
Overload the derivatives of the residuals with respect to a parameter to apply to the augmented syste...
AssemblyHandler * Assembly_handler_pt
The original assembly handler.
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
Vector< int > Count
A vector that is used to determine how many elements contribute to a particular equation. It is used to ensure that the global system is correctly formulated.
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
Problem * Problem_pt
Pointer to the problem.
cstr elem_len * i
Definition: cfortran.h:607
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
The linear-algebra-type solver does not make sense. The interface is deliberately broken...
virtual void get_inner_products(GeneralisedElement *const &elem_pt, Vector< std::pair< unsigned, unsigned > > const &history_index, Vector< double > &inner_product)
Compute the inner products of the given vector of pairs of history values over the element...
void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tra...
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Get the residuals.
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Get the global equation number of the local unknown.
int bifurcation_type() const
Indicate that we are tracking a fold bifurcation by returning 1.
BlockPitchForkLinearSolver(LinearSolver *const linear_solver_pt)
Constructor, inherits the original linear solver.
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Return the contribution to the residuals of the element elem_pt by using the derivatives.
The Problem class.
Definition: problem.h:152
Vector< int > Count
A vector that is used to determine how many elements contribute to a particular equation. It is used to ensure that the global system is correctly formulated.
A class that is used to assemble the residuals in parallel by overloading the get_all_vectors_and_mat...
Problem * Problem_pt
Pointer to the problem, used in the resolve.
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
Vector< double > Phi
A constant vector used to ensure that the null vector is not trivial.
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
virtual ~AssemblyHandler()
Empty virtual destructor.
virtual void get_all_vectors_and_matrices(GeneralisedElement *const &elem_pt, Vector< Vector< double > > &vec, Vector< DenseMatrix< double > > &matrix)
Calculate all desired vectors and matrices provided by the element elem_pt.
virtual void get_hessian_vector_products(GeneralisedElement *const &elem_pt, Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Calculate the product of the Hessian (derivative of Jacobian with respect to all variables) an eigenv...
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
The linear-algebra-type solver does not make sense.
double * Parameter_pt
Pointer to the parameter.
void synchronise()
Function that is used to perform any synchronisation required during the solution.
char t
Definition: cfortran.h:572
~HopfHandler()
Destructor return the problem to its original (non-augmented) state.
virtual void synchronise()
Function that is used to perform any synchronisation required during the solution.
DoubleVector * G_pt
Pointer to the storage for the vector g (0 to n-1)
DoubleVector * E_pt
Pointer to the storage for the vector e.
double * bifurcation_parameter_pt() const
Return a pointer to the bifurcation parameter in bifurcation tracking problems.
Problem * Problem_pt
Pointer to the problem, used in the resolve.
Vector< double > C
A constant vector used to ensure that the null vector is not trivial.
LinearAlgebraDistribution * Augmented_dof_distribution_pt
The augmented distribution.
~EigenProblemHandler()
Empty virtual destructor.
virtual void dof_vector(GeneralisedElement *const &elem_pt, const unsigned &t, Vector< double > &dof)
Return vector of dofs at time level t in the element elem_pt.
unsigned Solve_which_system
Integer flag to indicate which system should be assembled. There are three possibilities. The full augmented system (0), the non-augmented jacobian system (1), a system in which the jacobian is augmented by 1 row and column to ensure that it is non-singular (2). See the enum above.
double * Parameter_pt
The value of the parameter.
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Use underlying AssemblyHandler to return the contribution to the residuals of the element elem_pt...
unsigned Ndof
Store the number of degrees of freedom in the non-augmented problem.
void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable".
ParallelResidualsHandler(AssemblyHandler *const &assembly_handler_pt)
Constructor, set the original assembly handler.
DoubleVector * Alpha_pt
Pointer to the storage for the vector alpha.
AugmentedBlockPitchForkLinearSolver(LinearSolver *const linear_solver_pt)
Constructor, inherits the original linear solver.
LinearAlgebraDistribution * Dof_distribution_pt
Store the original dof distribution.
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Get the residuals.
void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tra...
Vector< double > Y
Storage for the null vector.
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
DoubleVectorWithHaloEntries C
A constant vector used to ensure that the null vector is not trivial.
unsigned ndof(GeneralisedElement *const &elem_pt)
Use underlying assembly handler to return the number of degrees of freedom in the element elem_pt...
virtual void dof_pt_vector(GeneralisedElement *const &elem_pt, Vector< double * > &dof_pt)
Return vector of pointers to dofs in the element elem_pt.
double Sigma_real
Storage for the real shift.
void solve_full_system()
Solve non-block system.
DoubleVector * B_pt
Pointer to the storage for the vector b.
virtual void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tra...
DoubleVector * C_pt
Pointer to the storage for the vector c.
unsigned Ndof
Store the number of degrees of freedom in the non-augmented problem.
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. Again deliberately br...
ParameterDerivativeHandler(AssemblyHandler *const &assembly_handler_pt, double *const &parameter_pt)
Store the original assembly handler and parameter.
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
unsigned ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
The linear-algebra-type solver does not make sense. The interface is deliberately broken...
int bifurcation_type() const
Indicate that we are tracking a Hopf bifurcation by returning 3.
DoubleVectorWithHaloEntries Y
Storage for the null vector.
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.
~AugmentedBlockPitchForkLinearSolver()
Destructor: clean up the allocated memory.
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
The linear-algebra-type solver does not make sense. The interface is deliberately broken...
void get_djacobian_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Overload the derivative of the residuals and jacobian with respect to a parameter so that it breaks...
bool Distributed
Boolean to indicate whether the problem is distributed.
virtual unsigned ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
void get_djacobian_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Overload the derivative of the residuals and jacobian with respect to a parameter so that it breaks...
void get_hessian_vector_products(GeneralisedElement *const &elem_pt, Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Overload the hessian vector product function so that it breaks.
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.Pass through to the orig...
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
~ParallelResidualsHandler()
Empty virtual destructor.
virtual void get_inner_product_vectors(GeneralisedElement *const &elem_pt, Vector< unsigned > const &history_index, Vector< Vector< double > > &inner_product_vector)
Compute the vectors that when taken as a dot product with other history values give the inner product...
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
unsigned ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
double * Parameter_pt
Storage for the pointer to the parameter.
~ExplicitTimeStepHandler()
Empty virtual destructor.
DoubleVector * E_pt
Pointer to the storage for the vector e.
virtual void get_djacobian_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Calculate the derivative of the residuals and jacobian with respect to a parameter.
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Use underlying AssemblyHandler to return the global equation number of the local unknown ieqn_local i...
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
FoldHandler(Problem *const &problem_pt, double *const &parameter_pt)
Constructor: initialise the fold handler, by setting initial guesses for Y, Phi and calculating count...
DoubleVector * D_pt
Pointer to the storage for the vector d.
void solve_for_two_rhs(Problem *const &problem_pt, DoubleVector &result, const DoubleVector &rhs2, DoubleVector &result2)
Solve for two right hand sides.
unsigned Solve_which_system
Integer flag to indicate which system should be assembled. There are three possibilities. The full augmented system (0), the non-augmented jacobian system (1), a system in which the jacobian is augmented by 1 row and column to ensure that it is non-singular (2). See the enum above.
DoubleVector * Alpha_pt
Pointer to the storage for the vector alpha.
AssemblyHandler()
Empty constructor.
AugmentedBlockFoldLinearSolver(LinearSolver *const linear_solver_pt)
Constructor, inherits the original linear solver.
void get_dresiduals_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam)
Overload the derivatives of the residuals with respect to a parameter to apply to the augmented syste...
void get_all_vectors_and_matrices(GeneralisedElement *const &elem_pt, Vector< Vector< double > > &vec, Vector< DenseMatrix< double > > &matrix)
Calculate all desired vectors and matrices provided by the element elem_pt.
A class that is used to define the functions used to assemble the elemental contributions to the mass...
virtual double & local_problem_dof(Problem *const &problem_pt, const unsigned &t, const unsigned &i)
Return the t-th level of storage associated with the i-th (local) dof stored in the problem...
void solve_block_system()
Set to solve the block system.
~FoldHandler()
Destructor, return the problem to its original state before the augmented system was added...
void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that.
unsigned ndof(GeneralisedElement *const &elem_pt)
Get the number of elemental degrees of freedom.
~AugmentedBlockFoldLinearSolver()
Destructor: clean up the allocated memory.
A class that is used to define the functions used to assemble the elemental contributions to the resi...
void solve_full_system()
Solve non-block system.
unsigned Solve_which_system
Integer flag to indicate which system should be assembled. There are three possibilities. The full augmented system (0), the non-augmented jacobian system (1), and complex system (2), where the matrix is a combination of the jacobian and mass matrices.
void get_djacobian_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Overload the derivative of the residuals and jacobian with respect to a parameter so that it breaks...
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Return the contribution to the residuals of the element elem_pt This is deliberately broken in our ei...
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
virtual void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Return the contribution to the residuals of the element elem_pt.
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
The linear-algebra-type solver does not make sense. The interface is deliberately broken...
HopfHandler(Problem *const &problem_pt, double *const &parameter_pt)
Constructor.
Problem * Problem_pt
Pointer to the problem, used in the resolve.
AssemblyHandler * Assembly_handler_pt
The original assembly handler.
virtual double * bifurcation_parameter_pt() const
Return a pointer to the bifurcation parameter in bifurcation tracking problems.
~PitchForkHandler()
Destructor return the problem to its original (non-augmented) state.
A class that is used to define the functions used when assembling the derivatives of the residuals wi...
void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable".
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.
void solve_complex_system()
Set to solve the complex system.
Problem * Problem_pt
Pointer to the problem.
unsigned ndof(GeneralisedElement *const &elem_pt)
virtual 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.
int bifurcation_type() const
Indicate that we are tracking a pitchfork bifurcation by returning 2.
Problem * Problem_pt
Pointer to the problem.
double * Parameter_pt
Storage for the pointer to the parameter.
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
The linear-algebra-type solver does not make sense. The interface is deliberately broken...
~BlockHopfLinearSolver()
Destructor: clean up the allocated memory.
unsigned ndof(GeneralisedElement *const &elem_pt)
Get the number of elemental degrees of freedom.
void solve_full_system()
Solve non-block system.
void solve_augmented_block_system()
Set to solve the augmented block system.
DoubleVectorWithHaloEntries Psi
A constant vector that is specifies the symmetry being broken.
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Get the global equation number of the local unknown.
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
virtual 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.
void get_hessian_vector_products(GeneralisedElement *const &elem_pt, Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Overload the hessian vector product function so that it breaks.
double * bifurcation_parameter_pt() const
Return a pointer to the bifurcation parameter in bifurcation tracking problems.
double * bifurcation_parameter_pt() const
Return a pointer to the bifurcation parameter in bifurcation tracking problems.
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
The linear-algebra-type solver does not make sense.
Vector< unsigned > Global_eqn_number
A vector that is used to map the global equations to their actual location in a distributed problem...
double Omega
The critical frequency of the bifurcation.
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
unsigned global_eqn_number(const unsigned &i)
Function that is used to return map the global equations using the simplistic numbering scheme into t...
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Return the contribution to the residuals of the element elem_pt This is deliberately broken in our ei...
unsigned ndof(GeneralisedElement *const &elem_pt)
Get the number of elemental degrees of freedom.
Vector< double > Phi
The real part of the null vector.
A class that is used to define the functions used to assemble and invert the mass matrix when taking ...
EigenProblemHandler(const double &sigma_real)
Constructor, sets the value of the real shift.
ExplicitTimeStepHandler()
Empty Constructor.
void solve_standard_system()
Set to solve the standard system.
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Get the global equation number of the local unknown.
void get_all_vectors_and_matrices(GeneralisedElement *const &elem_pt, Vector< Vector< double > > &vec, Vector< DenseMatrix< double > > &matrix)
Calculate all desired vectors and matrices provided by the element elem_pt This function calls only t...
void solve_augmented_block_system()
Set to solve the augmented block system.
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
virtual int bifurcation_type() const
Return an unsigned integer to indicate whether the handler is a bifurcation tracking handler...
Abstract base class for matrices of doubles – adds abstract interfaces for solving, LU decomposition and multiplication by vectors.
Definition: matrices.h:275
unsigned Ndof
Store the number of degrees of freedom in the non-augmented problem.
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Get the residuals.
void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable".
PitchForkHandler(Problem *const &problem_pt, AssemblyHandler *const &assembly_handler_pt, double *const &parameter_pt, const DoubleVector &symmetry_vector)
Constructor, initialise the systems.
const double & omega() const
Return the frequency of the bifurcation.
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. Again deliberately br...
Vector< double > Psi
The imaginary part of the null vector.
BlockHopfLinearSolver(LinearSolver *const linear_solver_pt)
Constructor, inherits the original linear solver.
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
The linear-algebra-type solver does not make sense. The interface is deliberately broken...
void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Use underlying AssemblyHandler to Calculate the elemental Jacobian matrix "d equation / d variable" ...
void get_hessian_vector_products(GeneralisedElement *const &elem_pt, Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Overload the hessian vector product function so that it breaks.
virtual void get_dresiduals_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam)
Calculate the derivative of the residuals with respect to a parameter.