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// 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// A class that is used to assemble the linear systems solved
27// by oomph-lib.
28
29// Include guards to prevent multiple inclusion of this header
30#ifndef OOMPH_ASSEMBLY_HANDLER_CLASS_HEADER
31#define OOMPH_ASSEMBLY_HANDLER_CLASS_HEADER
32
33// Config header generated by autoconfig
34#ifdef HAVE_CONFIG_H
35#include <oomph-lib-config.h>
36#endif
37
38// OOMPH-LIB headers
39#include "matrices.h"
40#include "linear_solver.h"
42
43namespace oomph
44{
45 // Forward class definition of the element
46 class GeneralisedElement;
47
48 // Forward class definition of the problem
49 class Problem;
50
51 //=============================================================
52 /// A class that is used to define the functions used to
53 /// assemble the elemental contributions to the
54 /// residuals vector and Jacobian matrix that define
55 /// the problem being solved.
56 /// The main use
57 /// of this class is to assemble and solve the augmented systems
58 /// used in bifurcation detection and tracking. The default
59 /// implementation merely calls the underlying elemental
60 /// functions with no augmentation.
61 //===============================================================
63 {
64 public:
65 /// Empty constructor
67
68 /// Return the number of degrees of freedom in the element elem_pt
69 virtual unsigned ndof(GeneralisedElement* const& elem_pt);
70
71 /// Return vector of dofs at time level t in the element elem_pt
72 virtual void dof_vector(GeneralisedElement* const& elem_pt,
73 const unsigned& t,
74 Vector<double>& dof);
75
76 /// Return vector of pointers to dofs in the element elem_pt
77 virtual void dof_pt_vector(GeneralisedElement* const& elem_pt,
78 Vector<double*>& dof_pt);
79
80 /// Return the t-th level of storage associated with the i-th
81 /// (local) dof stored in the problem
82 virtual double& local_problem_dof(Problem* const& problem_pt,
83 const unsigned& t,
84 const unsigned& i);
85
86 /// Return the global equation number of the local unknown ieqn_local
87 /// in elem_pt.
88 virtual unsigned long eqn_number(GeneralisedElement* const& elem_pt,
89 const unsigned& ieqn_local);
90
91 /// Return the contribution to the residuals of the element elem_pt
92 virtual void get_residuals(GeneralisedElement* const& elem_pt,
93 Vector<double>& residuals);
94
95 /// Calculate the elemental Jacobian matrix "d equation
96 /// / d variable" for elem_pt.
97 virtual void get_jacobian(GeneralisedElement* const& elem_pt,
98 Vector<double>& residuals,
99 DenseMatrix<double>& jacobian);
100
101 /// Calculate all desired vectors and matrices
102 /// provided by the element elem_pt.
103 virtual void get_all_vectors_and_matrices(
104 GeneralisedElement* const& elem_pt,
106 Vector<DenseMatrix<double>>& matrix);
107
108 /// Calculate the derivative of the residuals with respect to
109 /// a parameter
110 virtual void get_dresiduals_dparameter(GeneralisedElement* const& elem_pt,
111 double* const& parameter_pt,
112 Vector<double>& dres_dparam);
113
114 /// Calculate the derivative of the residuals and jacobian
115 /// with respect to a parameter
116 virtual void get_djacobian_dparameter(GeneralisedElement* const& elem_pt,
117 double* const& parameter_pt,
118 Vector<double>& dres_dparam,
119 DenseMatrix<double>& djac_dparam);
120
121 /// Calculate the product of the Hessian (derivative of Jacobian with
122 /// respect to all variables) an eigenvector, Y, and
123 /// other specified vectors, C
124 /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}
125 virtual void get_hessian_vector_products(GeneralisedElement* const& elem_pt,
126 Vector<double> const& Y,
127 DenseMatrix<double> const& C,
128 DenseMatrix<double>& product);
129
130
131 /// Return an unsigned integer to indicate whether the
132 /// handler is a bifurcation tracking handler. The default
133 /// is zero (not)
134 virtual int bifurcation_type() const
135 {
136 return 0;
137 }
138
139 /// Return a pointer to the
140 /// bifurcation parameter in bifurcation tracking problems
141 virtual double* bifurcation_parameter_pt() const;
142
143 /// 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 /// Compute the inner products of the given vector of pairs of
148 /// history values over the element.
149 virtual void get_inner_products(
150 GeneralisedElement* const& elem_pt,
151 Vector<std::pair<unsigned, unsigned>> const& history_index,
152 Vector<double>& inner_product);
153
154 /// 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(
157 GeneralisedElement* const& elem_pt,
158 Vector<unsigned> const& history_index,
159 Vector<Vector<double>>& inner_product_vector);
160
161#ifdef OOMPH_HAS_MPI
162
163 /// Function that is used to perform any synchronisation
164 /// required during the solution
165 virtual void synchronise() {}
166#endif
167
168 /// Empty virtual destructor
169 virtual ~AssemblyHandler() {}
170 };
171
172
173 //=============================================================
174 /// A class that is used to define the functions used to
175 /// assemble and invert the mass matrix when taking an explicit
176 /// timestep. The idea is simply to replace the jacobian matrix
177 /// with the mass matrix and then our standard linear solvers
178 /// will solve the required system
179 //===============================================================
181 {
182 public:
183 /// Empty Constructor
185
186 /// Return the number of degrees of freedom in the element elem_pt
187 unsigned ndof(GeneralisedElement* const& elem_pt);
188
189 /// Return the global equation number of the local unknown ieqn_local
190 /// in elem_pt.
191 unsigned long eqn_number(GeneralisedElement* const& elem_pt,
192 const unsigned& ieqn_local);
193
194 /// Return the contribution to the residuals of the element elem_pt
195 /// This is deliberately broken in our eigenproblem
196 void get_residuals(GeneralisedElement* const& elem_pt,
197 Vector<double>& residuals);
198
199 /// Calculate the elemental Jacobian matrix "d equation
200 /// / d variable" for elem_pt. Again deliberately broken in the eigenproblem
201 void get_jacobian(GeneralisedElement* const& elem_pt,
202 Vector<double>& residuals,
203 DenseMatrix<double>& jacobian);
204
205 /// Calculate all desired vectors and matrices
206 /// provided by the element elem_pt.
209 Vector<DenseMatrix<double>>& matrix);
210
211 /// Empty virtual destructor
213 };
214
215
216 //=============================================================
217 /// A class that is used to define the functions used to
218 /// assemble the elemental contributions to the
219 /// mass matrix and jacobian (stiffness) matrix that define
220 /// a generalised eigenproblem.
221 //===============================================================
223 {
224 /// Storage for the real shift
226
227 public:
228 /// Constructor, sets the value of the real shift
229 EigenProblemHandler(const double& sigma_real) : Sigma_real(sigma_real) {}
230
231 /// Return the number of degrees of freedom in the element elem_pt
232 unsigned ndof(GeneralisedElement* const& elem_pt);
233
234 /// Return the global equation number of the local unknown ieqn_local
235 /// in elem_pt.
236 unsigned long eqn_number(GeneralisedElement* const& elem_pt,
237 const unsigned& ieqn_local);
238
239 /// Return the contribution to the residuals of the element elem_pt
240 /// This is deliberately broken in our eigenproblem
241 void get_residuals(GeneralisedElement* const& elem_pt,
242 Vector<double>& residuals);
243
244 /// Calculate the elemental Jacobian matrix "d equation
245 /// / d variable" for elem_pt. Again deliberately broken in the eigenproblem
246 void get_jacobian(GeneralisedElement* const& elem_pt,
247 Vector<double>& residuals,
248 DenseMatrix<double>& jacobian);
249
250 /// Calculate all desired vectors and matrices
251 /// provided by the element elem_pt.
254 Vector<DenseMatrix<double>>& matrix);
255
256 /// Empty virtual destructor
258 };
259
260
261 //=============================================================
262 /// A class that is used to assemble the residuals in
263 /// parallel by overloading the get_all_vectors_and_matrices,
264 /// so that only the residuals are returned. This ensures that
265 /// the (moderately complex) distributed parallel assembly
266 /// loops are only in one place.
267 //===============================================================
269 {
270 /// The original assembly handler
272
273 public:
274 /// Constructor, set the original assembly handler
275 ParallelResidualsHandler(AssemblyHandler* const& assembly_handler_pt)
276 : Assembly_handler_pt(assembly_handler_pt)
277 {
278 }
279
280 /// Use underlying assembly handler to return the number of
281 /// degrees of freedom in the element elem_pt
282 unsigned ndof(GeneralisedElement* const& elem_pt)
283 {
284 return Assembly_handler_pt->ndof(elem_pt);
285 }
286
287 /// Use underlying AssemblyHandler to return the
288 /// global equation number of the local unknown ieqn_local in elem_pt.
289 unsigned long eqn_number(GeneralisedElement* const& elem_pt,
290 const unsigned& ieqn_local)
291 {
292 return Assembly_handler_pt->eqn_number(elem_pt, ieqn_local);
293 }
294
295 /// Use underlying AssemblyHandler to return
296 /// the contribution to the residuals of the element elem_pt
297 void get_residuals(GeneralisedElement* const& elem_pt,
298 Vector<double>& residuals)
299 {
300 Assembly_handler_pt->get_residuals(elem_pt, residuals);
301 }
302
303
304 /// Use underlying AssemblyHandler to
305 /// Calculate the elemental Jacobian matrix "d equation
306 /// / d variable" for elem_pt.
307 void get_jacobian(GeneralisedElement* const& elem_pt,
308 Vector<double>& residuals,
309 DenseMatrix<double>& jacobian)
310 {
311 Assembly_handler_pt->get_jacobian(elem_pt, residuals, jacobian);
312 }
313
314
315 /// Calculate all desired vectors and matrices
316 /// provided by the element elem_pt
317 /// This function calls only the get_residuals function associated
318 /// with the original assembly handler
322 {
323 Assembly_handler_pt->get_residuals(elem_pt, vec[0]);
324 }
325
326 /// Empty virtual destructor
328 };
329
330
331 //==========================================================================
332 /// A class that is used to define the functions used when assembling
333 /// the derivatives of the residuals with respect to a parameter.
334 /// The idea is to replace get_residuals with get_dresiduals_dparameter with
335 /// a particular parameter and assembly handler that are passed on
336 /// assembly.
337 //==========================================================================
339 {
340 /// The value of the parameter
342
343 /// The original assembly handler
345
346 public:
347 /// Store the original assembly handler and parameter
348 ParameterDerivativeHandler(AssemblyHandler* const& assembly_handler_pt,
349 double* const& parameter_pt)
350 : Parameter_pt(parameter_pt), Assembly_handler_pt(assembly_handler_pt)
351 {
352 }
353
354 /// Return the number of degrees of freedom in the element elem_pt
355 /// Pass through to the original assembly handler
356 unsigned ndof(GeneralisedElement* const& elem_pt)
357 {
358 return Assembly_handler_pt->ndof(elem_pt);
359 }
360
361 /// Return the global equation number of the local unknown ieqn_local
362 /// in elem_pt.Pass through to the original assembly handler
363 unsigned long eqn_number(GeneralisedElement* const& elem_pt,
364 const unsigned& ieqn_local)
365 {
366 return Assembly_handler_pt->eqn_number(elem_pt, ieqn_local);
367 }
368
369 /// Return the contribution to the residuals of the element elem_pt
370 /// by using the derivatives
371 void get_residuals(GeneralisedElement* const& elem_pt,
372 Vector<double>& residuals)
373 {
375 elem_pt, Parameter_pt, residuals);
376 }
377
378
379 /// Calculate the elemental Jacobian matrix "d equation
380 /// / d variable" for elem_pt.
381 /// Overloaded to return the derivatives wrt the parameter
382 void get_jacobian(GeneralisedElement* const& elem_pt,
383 Vector<double>& residuals,
384 DenseMatrix<double>& jacobian)
385 {
387 elem_pt, Parameter_pt, residuals, jacobian);
388 }
389 };
390
391
392 //========================================================================
393 /// A custom linear solver class that is used to solve a block-factorised
394 /// version of the Fold bifurcation detection problem.
395 //========================================================================
397 {
398 /// Pointer to the original linear solver
400
401 // Pointer to the problem, used in the resolve
403
404 /// Pointer to the storage for the vector alpha
406
407 /// Pointer to the storage for the vector e
409
410 public:
411 /// Constructor, inherits the original linear solver
414 {
415 }
416
417 /// Destructor: clean up the allocated memory
419
420 /// The solve function uses the block factorisation
421 void solve(Problem* const& problem_pt, DoubleVector& result);
422
423 /// The linear-algebra-type solver does not make sense.
424 /// The interface is deliberately broken
425 void solve(DoubleMatrixBase* const& matrix_pt,
426 const DoubleVector& rhs,
427 DoubleVector& result)
428 {
429 throw OomphLibError(
430 "Linear-algebra interface does not make sense for this linear solver\n",
431 OOMPH_CURRENT_FUNCTION,
432 OOMPH_EXCEPTION_LOCATION);
433 }
434
435 /// The linear-algebra-type solver does not make sense.
436 /// The interface is deliberately broken
437 void solve(DoubleMatrixBase* const& matrix_pt,
438 const Vector<double>& rhs,
439 Vector<double>& result)
440 {
441 throw OomphLibError(
442 "Linear-algebra interface does not make sense for this linear solver\n",
443 OOMPH_CURRENT_FUNCTION,
444 OOMPH_EXCEPTION_LOCATION);
445 }
446
447 /// The resolve function also uses the block factorisation
448 void resolve(const DoubleVector& rhs, DoubleVector& result);
449
450 /// Access function to the original linear solver
452 {
453 return Linear_solver_pt;
454 }
455 };
456
457
458 //===================================================================
459 /// A class that is used to assemble the augmented system that defines
460 /// a fold (saddle-node) or limit point. The "standard" problem must
461 /// be a function of a global paramter \f$\lambda\f$, and a
462 /// solution is \f$R(u,\lambda) = 0 \f$ , where \f$ u \f$ are the unknowns
463 /// in the problem. A limit point is formally specified by the augmented
464 /// system of size \f$ 2N+1 \f$
465 /// \f[ R(u,\lambda) = 0, \f]
466 /// \f[ Jy = 0, \f]
467 /// \f[\phi\cdot y = 1. \f]
468 /// In the above \f$ J \f$ is the usual Jacobian matrix, \f$ dR/du \f$, and
469 /// \f$ \phi \f$ is a constant vector that is chosen to
470 /// ensure that the null vector, \f$ y \f$, is not trivial.
471 //======================================================================
473 {
475
476 /// A little private enum to determine whether we are solving
477 /// the block system or not
478 enum
479 {
483 };
484
485 /// Integer flag to indicate which system should be assembled.
486 /// There are three possibilities. The full augmented system (0),
487 /// the non-augmented jacobian system (1), a system in which the
488 /// jacobian is augmented by 1 row and column to ensure that it is
489 /// non-singular (2). See the enum above
491
492 /// Pointer to the problem
494
495 /// Store the number of degrees of freedom in the non-augmented
496 /// problem
497 unsigned Ndof;
498
499 /// A constant vector used to ensure that the null vector
500 /// is not trivial
502
503 /// Storage for the null vector
505
506 /// A vector that is used to determine how many elements
507 /// contribute to a particular equation. It is used to ensure
508 /// that the global system is correctly formulated.
510
511 /// Storage for the pointer to the parameter
513
514 public:
515 /// Constructor:
516 /// initialise the fold handler, by setting initial guesses
517 /// for Y, Phi and calculating count. If the system changes, a new
518 /// fold handler must be constructed
519 FoldHandler(Problem* const& problem_pt, double* const& parameter_pt);
520
521
522 /// Constructor in which initial eigenvector can be passed
523 FoldHandler(Problem* const& problem_pt,
524 double* const& parameter_pt,
525 const DoubleVector& eigenvector);
526
527 /// Constructor in which initial eigenvector
528 /// and normalisation can be passed
529 FoldHandler(Problem* const& problem_pt,
530 double* const& parameter_pt,
531 const DoubleVector& eigenvector,
532 const DoubleVector& normalisation);
533
534
535 /// Destructor, return the problem to its original state
536 /// before the augmented system was added
537 ~FoldHandler();
538
539 /// Get the number of elemental degrees of freedom
540 unsigned ndof(GeneralisedElement* const& elem_pt);
541
542 /// Get the global equation number of the local unknown
543 unsigned long eqn_number(GeneralisedElement* const& elem_pt,
544 const unsigned& ieqn_local);
545
546 /// Get the residuals
547 void get_residuals(GeneralisedElement* const& elem_pt,
548 Vector<double>& residuals);
549
550 /// Calculate the elemental Jacobian matrix "d equation
551 /// / d variable".
552 void get_jacobian(GeneralisedElement* const& elem_pt,
553 Vector<double>& residuals,
554 DenseMatrix<double>& jacobian);
555
556 /// Overload the derivatives of the residuals with respect to
557 /// a parameter to apply to the augmented system
559 double* const& parameter_pt,
560 Vector<double>& dres_dparam);
561
562 /// Overload the derivative of the residuals and jacobian
563 /// with respect to a parameter so that it breaks
564 void get_djacobian_dparameter(GeneralisedElement* const& elem_pt,
565 double* const& parameter_pt,
566 Vector<double>& dres_dparam,
567 DenseMatrix<double>& djac_dparam);
568
569 /// Overload the hessian vector product function so that
570 /// it breaks
572 Vector<double> const& Y,
573 DenseMatrix<double> const& C,
574 DenseMatrix<double>& product);
575
576 /// Indicate that we are tracking a fold bifurcation by returning 1
578 {
579 return 1;
580 }
581
582 /// Return a pointer to the
583 /// bifurcation parameter in bifurcation tracking problems
585 {
586 return Parameter_pt;
587 }
588
589 /// Return the eigenfunction(s) associated with the bifurcation that
590 /// has been detected in bifurcation tracking problems
591 void get_eigenfunction(Vector<DoubleVector>& eigenfunction);
592
593 /// Set to solve the augmented block system
595
596 /// Set to solve the block system
597 void solve_block_system();
598
599 /// Solve non-block system
600 void solve_full_system();
601 };
602
603
604 //========================================================================
605 /// A custom linear solver class that is used to solve a block-factorised
606 /// version of the PitchFork bifurcation detection problem.
607 //========================================================================
609 {
610 /// Pointer to the original linear solver
612
613 /// Pointer to the problem, used in the resolve
615
616 /// Pointer to the storage for the vector b
618
619 /// Pointer to the storage for the vector c
621
622 /// Pointer to the storage for the vector d
624
625 /// Pointer to the storage for the vector of derivatives with respect
626 /// to the bifurcation parameter
628
629 public:
630 /// Constructor, inherits the original linear solver
633 Problem_pt(0),
634 B_pt(0),
635 C_pt(0),
636 D_pt(0),
638 {
639 }
640
641 /// Destructor: clean up the allocated memory
643
644 /// The solve function uses the block factorisation
645 void solve(Problem* const& problem_pt, DoubleVector& result);
646
647 /// The linear-algebra-type solver does not make sense.
648 /// The interface is deliberately broken
649 void solve(DoubleMatrixBase* const& matrix_pt,
650 const DoubleVector& rhs,
651 DoubleVector& result)
652 {
653 throw OomphLibError(
654 "Linear-algebra interface does not make sense for this linear solver\n",
655 OOMPH_CURRENT_FUNCTION,
656 OOMPH_EXCEPTION_LOCATION);
657 }
658
659 /// The linear-algebra-type solver does not make sense.
660 /// The interface is deliberately broken
661 void solve(DoubleMatrixBase* const& matrix_pt,
662 const Vector<double>& rhs,
663 Vector<double>& result)
664 {
665 throw OomphLibError(
666 "Linear-algebra interface does not make sense for this linear solver\n",
667 OOMPH_CURRENT_FUNCTION,
668 OOMPH_EXCEPTION_LOCATION);
669 }
670
671
672 /// The resolve function also uses the block factorisation
673 void resolve(const DoubleVector& rhs, DoubleVector& result);
674
675 /// Access function to the original linear solver
677 {
678 return Linear_solver_pt;
679 }
680 };
681
682
683 //========================================================================
684 /// A custom linear solver class that is used to solve a block-factorised
685 /// version of the PitchFork bifurcation detection problem.
686 //========================================================================
688 {
689 /// Pointer to the original linear solver
691
692 /// Pointer to the problem, used in the resolve
694
695 /// Pointer to the storage for the vector alpha
697
698 /// Pointer to the storage for the vector e
700
701 public:
702 /// Constructor, inherits the original linear solver
705 {
706 }
707
708 /// Destructor: clean up the allocated memory
710
711 /// The solve function uses the block factorisation
712 void solve(Problem* const& problem_pt, DoubleVector& result);
713
714 /// The linear-algebra-type solver does not make sense.
715 /// The interface is deliberately broken
716 void solve(DoubleMatrixBase* const& matrix_pt,
717 const DoubleVector& rhs,
718 DoubleVector& result)
719 {
720 throw OomphLibError(
721 "Linear-algebra interface does not make sense for this linear solver\n",
722 OOMPH_CURRENT_FUNCTION,
723 OOMPH_EXCEPTION_LOCATION);
724 }
725
726 /// The linear-algebra-type solver does not make sense.
727 /// The interface is deliberately broken
728 void solve(DoubleMatrixBase* const& matrix_pt,
729 const Vector<double>& rhs,
730 Vector<double>& result)
731 {
732 throw OomphLibError(
733 "Linear-algebra interface does not make sense for this linear solver\n",
734 OOMPH_CURRENT_FUNCTION,
735 OOMPH_EXCEPTION_LOCATION);
736 }
737
738
739 /// The resolve function also uses the block factorisation
740 void resolve(const DoubleVector& rhs, DoubleVector& result);
741
742 /// Access function to the original linear solver
744 {
745 return Linear_solver_pt;
746 }
747 };
748
749 //========================================================================
750 /// A class that is used to assemble the augmented system that defines
751 /// a pitchfork (symmetry-breaking) bifurcation. The "standard" problem
752 /// must be a function of a global parameter \f$ \lambda \f$ and a solution
753 /// is \f$R(u,\lambda) = 0\f$, where \f$u\f$ are the unknowns in the
754 /// problem. A pitchfork bifurcation may be specified by the augmented
755 /// system of size \f$2N+2\f$.
756 /// \f[ R(u,\lambda) + \sigma \psi = 0,\f]
757 /// \f[ Jy = 0,\f]
758 /// \f[ \langle u, \phi \rangle = 0, \f]
759 /// \f[ \phi \cdot y = 1.\f]
760 /// In the abovem \f$J\f$ is the usual Jacobian matrix, \f$dR/du\f$
761 /// and \f$\phi\f$ is a constant vector that is chosen to ensure that
762 /// the null vector, \f$y\f$, is not trivial.
763 /// Here \f$\sigma \f$ is a slack variable that is used to enforce the
764 /// constraint \f$ \langle u, \phi \rangle = 0 \f$ --- the inner product
765 /// of the solution with the chosen symmetry vector is zero. At the
766 /// bifurcation \f$ \sigma \f$ should be very close to zero. Note that
767 /// this formulation means that any odd symmetry in the problem must
768 /// be about zero. Moreover, we use the dot product of two vectors to
769 /// calculate the inner product, rather than an integral over the
770 /// domain and so the mesh must be symmetric in the appropriate directions.
771 //========================================================================
773 {
776
777 /// A little private enum to determine whether we are solving
778 /// the block system or not
779 enum
780 {
784 };
785
786 /// Integer flag to indicate which system should be assembled.
787 /// There are three possibilities. The full augmented system (0),
788 /// the non-augmented jacobian system (1), a system in which the
789 /// jacobian is augmented by 1 row and column to ensure that it is
790 /// non-singular (2). See the enum above
792
793 /// Pointer to the problem
795
796 /// Pointer to the underlying (original) assembly handler
798
799 /// Store the number of degrees of freedom in the non-augmented
800 /// problem
801 unsigned Ndof;
802
803 /// Store the original dof distribution
805
806 /// The augmented distribution
808
809 /// A slack variable used to specify the amount of antisymmetry
810 /// in the solution.
811 double Sigma;
812
813 /// Storage for the null vector
815
816 /// A constant vector that is specifies the symmetry being broken
818
819 /// A constant vector used to ensure that the null vector
820 /// is not trivial
822
823 /// A vector that is used to determine how many elements
824 /// contribute to a particular equation. It is used to ensure
825 /// that the global system is correctly formulated.
826 /// This should really be an integer, but its double so that
827 /// the distribution can be used
829
830 /// A vector that is used to map the global equations to their
831 /// actual location in a distributed problem
833
834 /// Storage for the pointer to the parameter
836
837 // The total number of elements in the problem
838 unsigned Nelement;
839
840#ifdef OOMPH_HAS_MPI
841 /// Boolean to indicate whether the problem is distributed
843#endif
844
845 /// Function that is used to return map the global equations
846 /// using the simplistic numbering scheme into the actual distributed
847 /// scheme
848 inline unsigned global_eqn_number(const unsigned& i)
849 {
850#ifdef OOMPH_HAS_MPI
851 // If the problem is distributed I have to do something
852 if (Distributed)
853 {
854 return Global_eqn_number[i];
855 }
856 // Otherwise it's just i
857 else
858#endif
859 {
860 return i;
861 }
862 }
863
864 public:
865 /// Constructor, initialise the systems
866 PitchForkHandler(Problem* const& problem_pt,
867 AssemblyHandler* const& assembly_handler_pt,
868 double* const& parameter_pt,
869 const DoubleVector& symmetry_vector);
870
871 /// Destructor, return the problem to its original state,
872 /// before the augmented system was added
874
875 // Has this been called
876
877 /// Get the number of elemental degrees of freedom
878 unsigned ndof(GeneralisedElement* const& elem_pt);
879
880 /// Get the global equation number of the local unknown
881 unsigned long eqn_number(GeneralisedElement* const& elem_pt,
882 const unsigned& ieqn_local);
883
884 /// Get the residuals
885 void get_residuals(GeneralisedElement* const& elem_pt,
886 Vector<double>& residuals);
887
888 /// Calculate the elemental Jacobian matrix "d equation
889 /// / d variable".
890 void get_jacobian(GeneralisedElement* const& elem_pt,
891 Vector<double>& residuals,
892 DenseMatrix<double>& jacobian);
893
894 /// Overload the derivatives of the residuals with respect to
895 /// a parameter to apply to the augmented system
897 double* const& parameter_pt,
898 Vector<double>& dres_dparam);
899
900 /// Overload the derivative of the residuals and jacobian
901 /// with respect to a parameter so that it breaks
902 void get_djacobian_dparameter(GeneralisedElement* const& elem_pt,
903 double* const& parameter_pt,
904 Vector<double>& dres_dparam,
905 DenseMatrix<double>& djac_dparam);
906
907 /// Overload the hessian vector product function so that
908 /// it breaks
910 Vector<double> const& Y,
911 DenseMatrix<double> const& C,
912 DenseMatrix<double>& product);
913
914
915 /// Indicate that we are tracking a pitchfork
916 /// bifurcation by returning 2
918 {
919 return 2;
920 }
921
922 /// Return a pointer to the
923 /// bifurcation parameter in bifurcation tracking problems
925 {
926 return Parameter_pt;
927 }
928
929 /// Return the eigenfunction(s) associated with the bifurcation that
930 /// has been detected in bifurcation tracking problems
931 void get_eigenfunction(Vector<DoubleVector>& eigenfunction);
932
933#ifdef OOMPH_HAS_MPI
934 /// Function that is used to perform any synchronisation
935 /// required during the solution
936 void synchronise();
937#endif
938
939
940 /// Set to solve the augmented block system
942
943 /// Set to solve the block system
944 void solve_block_system();
945
946 /// Solve non-block system
947 void solve_full_system();
948 };
949
950 //========================================================================
951 /// A custom linear solver class that is used to solve a block-factorised
952 /// version of the Hopf bifurcation detection problem.
953 //========================================================================
955 {
956 /// Pointer to the original linear solver
958
959 /// Pointer to the problem, used in the resolve
961
962 /// Pointer to the storage for the vector a
964
965 /// Pointer to the storage for the vector e (0 to n-1)
967
968 /// Pointer to the storage for the vector g (0 to n-1)
970
971 public:
972 /// Constructor, inherits the original linear solver
975 Problem_pt(0),
976 A_pt(0),
977 E_pt(0),
978 G_pt(0)
979 {
980 }
981
982 /// Destructor: clean up the allocated memory
984
985 /// Solve for two right hand sides
986 void solve_for_two_rhs(Problem* const& problem_pt,
987 DoubleVector& result,
988 const DoubleVector& rhs2,
989 DoubleVector& result2);
990
991 /// The solve function uses the block factorisation
992 void solve(Problem* const& problem_pt, DoubleVector& result);
993
994 /// The linear-algebra-type solver does not make sense.
995 /// The interface is deliberately broken
996 void solve(DoubleMatrixBase* const& matrix_pt,
997 const DoubleVector& rhs,
998 DoubleVector& result)
999 {
1000 throw OomphLibError(
1001 "Linear-algebra interface does not make sense for this linear solver\n",
1002 OOMPH_CURRENT_FUNCTION,
1003 OOMPH_EXCEPTION_LOCATION);
1004 }
1005
1006 /// The linear-algebra-type solver does not make sense.
1007 /// The interface is deliberately broken
1008 void solve(DoubleMatrixBase* const& matrix_pt,
1009 const Vector<double>& rhs,
1010 Vector<double>& result)
1011 {
1012 throw OomphLibError(
1013 "Linear-algebra interface does not make sense for this linear solver\n",
1014 OOMPH_CURRENT_FUNCTION,
1015 OOMPH_EXCEPTION_LOCATION);
1016 }
1017
1018
1019 /// The resolve function also uses the block factorisation
1020 void resolve(const DoubleVector& rhs, DoubleVector& result);
1021
1022 /// Access function to the original linear solver
1024 {
1025 return Linear_solver_pt;
1026 }
1027 };
1028
1029
1030 //===============================================================
1031 /// A class that is used to assemble the augmented system that defines
1032 /// a Hopf bifurcation. The "standard" problem
1033 /// must be a function of a global parameter \f$ \lambda \f$ and a solution
1034 /// is \f$ R(u,\lambda) = 0 \f$, where \f$ u \f$ are the unknowns in the
1035 /// problem. A Hopf bifurcation may be specified by the augmented
1036 /// system of size \f$ 3N+2 \f$.
1037 /// \f[ R(u,\lambda) = 0, \f]
1038 /// \f[ J\phi + \omega M \psi = 0, \f]
1039 /// \f[ J\psi - \omega M \phi = 0, \f]
1040 /// \f[ c \cdot \phi = 1. \f]
1041 /// \f[ c \cdot \psi = 0. \f]
1042 /// In the above \f$ J \f$ is the usual Jacobian matrix, \f$ dR/du \f$
1043 /// and \f$ M \f$ is the mass matrix that multiplies the time derivative
1044 /// terms. \f$ \phi + i\psi \f$ is the (complex) null vector of the complex
1045 /// matrix \f$ J - i\omega M \f$, where \f$ \omega \f$ is the critical
1046 /// frequency. \f$ c \f$ is a constant vector that is used to ensure that the
1047 /// null vector is non-trivial.
1048 //===========================================================================
1050 {
1052
1053 /// Integer flag to indicate which system should be assembled.
1054 /// There are three possibilities. The full augmented system (0),
1055 /// the non-augmented jacobian system (1), and complex
1056 /// system (2), where the matrix is a combination of the jacobian
1057 /// and mass matrices.
1059
1060 /// Pointer to the problem
1062
1063 /// Pointer to the parameter
1065
1066 /// Store the number of degrees of freedom in the non-augmented
1067 /// problem
1068 unsigned Ndof;
1069
1070 /// The critical frequency of the bifurcation
1071 double Omega;
1072
1073 /// The real part of the null vector
1075
1076 /// The imaginary part of the null vector
1078
1079 /// A constant vector used to ensure that the null vector is
1080 /// not trivial
1082
1083 /// A vector that is used to determine how many elements
1084 /// contribute to a particular equation. It is used to ensure
1085 /// that the global system is correctly formulated.
1087
1088 public:
1089 /// Constructor
1090 HopfHandler(Problem* const& problem_pt, double* const& parameter_pt);
1091
1092 /// Constructor with initial guesses for the frequency and null
1093 /// vectors, such as might be provided by an eigensolver
1094 HopfHandler(Problem* const& problem_pt,
1095 double* const& paramter_pt,
1096 const double& omega,
1097 const DoubleVector& phi,
1098 const DoubleVector& psi);
1099
1100 /// Destructor, return the problem to its original state,
1101 /// before the augmented system was added
1102 ~HopfHandler();
1103
1104 /// Get the number of elemental degrees of freedom
1105 unsigned ndof(GeneralisedElement* const& elem_pt);
1106
1107 /// Get the global equation number of the local unknown
1108 unsigned long eqn_number(GeneralisedElement* const& elem_pt,
1109 const unsigned& ieqn_local);
1110
1111 /// Get the residuals
1112 void get_residuals(GeneralisedElement* const& elem_pt,
1113 Vector<double>& residuals);
1114
1115 /// Calculate the elemental Jacobian matrix "d equation
1116 /// / d variable".
1117 void get_jacobian(GeneralisedElement* const& elem_pt,
1118 Vector<double>& residuals,
1119 DenseMatrix<double>& jacobian);
1120
1121 /// Overload the derivatives of the residuals with respect to
1122 /// a parameter to apply to the augmented system
1123 void get_dresiduals_dparameter(GeneralisedElement* const& elem_pt,
1124 double* const& parameter_pt,
1125 Vector<double>& dres_dparam);
1126
1127 /// Overload the derivative of the residuals and jacobian
1128 /// with respect to a parameter so that it breaks
1129 void get_djacobian_dparameter(GeneralisedElement* const& elem_pt,
1130 double* const& parameter_pt,
1131 Vector<double>& dres_dparam,
1132 DenseMatrix<double>& djac_dparam);
1133
1134 /// Overload the hessian vector product function so that
1135 /// it breaks
1137 Vector<double> const& Y,
1138 DenseMatrix<double> const& C,
1139 DenseMatrix<double>& product);
1140
1141 /// Indicate that we are tracking a Hopf
1142 /// bifurcation by returning 3
1144 {
1145 return 3;
1146 }
1147
1148 /// Return a pointer to the
1149 /// bifurcation parameter in bifurcation tracking problems
1151 {
1152 return Parameter_pt;
1153 }
1154
1155 /// Return the eigenfunction(s) associated with the bifurcation that
1156 /// has been detected in bifurcation tracking problems
1157 void get_eigenfunction(Vector<DoubleVector>& eigenfunction);
1158
1159 /// Return the frequency of the bifurcation
1160 const double& omega() const
1161 {
1162 return Omega;
1163 }
1164
1165 /// Set to solve the standard system
1166 void solve_standard_system();
1167
1168 /// Set to solve the complex system
1169 void solve_complex_system();
1170
1171 /// Solve non-block system
1172 void solve_full_system();
1173 };
1174
1175
1176} // namespace oomph
1177#endif
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...
virtual unsigned ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
virtual void synchronise()
Function that is used to perform any synchronisation required during the solution.
virtual int bifurcation_type() const
Return an unsigned integer to indicate whether the handler is a bifurcation tracking handler....
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...
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...
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.
virtual double * bifurcation_parameter_pt() const
Return a pointer to the bifurcation parameter in bifurcation tracking problems.
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.
virtual void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tra...
virtual void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Return the contribution to the residuals of the element elem_pt.
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.
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 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.
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.
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.
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.
AssemblyHandler()
Empty constructor.
virtual void dof_pt_vector(GeneralisedElement *const &elem_pt, Vector< double * > &dof_pt)
Return vector of pointers to dofs in the element elem_pt.
A custom linear solver class that is used to solve a block-factorised version of the Fold bifurcation...
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
AugmentedBlockFoldLinearSolver(LinearSolver *const linear_solver_pt)
Constructor, inherits the original linear solver.
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
DoubleVector * E_pt
Pointer to the storage for the vector e.
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
~AugmentedBlockFoldLinearSolver()
Destructor: clean up the allocated memory.
DoubleVector * Alpha_pt
Pointer to the storage for the vector alpha.
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 solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
The linear-algebra-type solver does not make sense. The interface is deliberately broken.
A custom linear solver class that is used to solve a block-factorised version of the PitchFork bifurc...
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
DoubleVector * Alpha_pt
Pointer to the storage for the vector alpha.
DoubleVector * E_pt
Pointer to the storage for the vector e.
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.
AugmentedBlockPitchForkLinearSolver(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.
~AugmentedBlockPitchForkLinearSolver()
Destructor: clean up the allocated memory.
Problem * Problem_pt
Pointer to the problem, used in the resolve.
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
A custom linear solver class that is used to solve a block-factorised version of the Hopf bifurcation...
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
DoubleVector * A_pt
Pointer to the storage for the vector a.
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 solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
Problem * Problem_pt
Pointer to the problem, used in the resolve.
~BlockHopfLinearSolver()
Destructor: clean up the allocated memory.
void solve_for_two_rhs(Problem *const &problem_pt, DoubleVector &result, const DoubleVector &rhs2, DoubleVector &result2)
Solve for two right hand sides.
DoubleVector * E_pt
Pointer to the storage for the vector e (0 to n-1)
BlockHopfLinearSolver(LinearSolver *const linear_solver_pt)
Constructor, inherits the original linear solver.
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.
DoubleVector * G_pt
Pointer to the storage for the vector g (0 to n-1)
A custom linear solver class that is used to solve a block-factorised version of the PitchFork bifurc...
DoubleVector * D_pt
Pointer to the storage for the vector d.
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
DoubleVector * dJy_dparam_pt
Pointer to the storage for the vector of derivatives with respect to the bifurcation parameter.
~BlockPitchForkLinearSolver()
Destructor: clean up the allocated memory.
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.
DoubleVector * B_pt
Pointer to the storage for the vector b.
Problem * Problem_pt
Pointer to the problem, used in the resolve.
BlockPitchForkLinearSolver(LinearSolver *const linear_solver_pt)
Constructor, inherits the original linear solver.
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
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.
DoubleVector * C_pt
Pointer to the storage for the vector c.
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
Abstract base class for matrices of doubles – adds abstract interfaces for solving,...
Definition: matrices.h:261
===================================================================== An extension of DoubleVector th...
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
A class that is used to define the functions used to assemble the elemental contributions to the mass...
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...
double Sigma_real
Storage for the real shift.
~EigenProblemHandler()
Empty virtual destructor.
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. Again deliberately bro...
EigenProblemHandler(const double &sigma_real)
Constructor, sets the value of the real shift.
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 ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
A class that is used to define the functions used to assemble and invert the mass matrix when taking ...
~ExplicitTimeStepHandler()
Empty virtual destructor.
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 ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
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.
ExplicitTimeStepHandler()
Empty Constructor.
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 bro...
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...
A class that is used to assemble the augmented system that defines a fold (saddle-node) or limit poin...
double * bifurcation_parameter_pt() const
Return a pointer to the bifurcation parameter in bifurcation tracking problems.
Problem * Problem_pt
Pointer to the problem.
Vector< int > Count
A vector that is used to determine how many elements contribute to a particular equation....
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Get the global equation number of the local unknown.
void solve_full_system()
Solve non-block system.
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.
unsigned Solve_which_system
Integer flag to indicate which system should be assembled. There are three possibilities....
void solve_augmented_block_system()
Set to solve the augmented block system.
Vector< double > Phi
A constant vector used to ensure that the null vector is not trivial.
unsigned Ndof
Store the number of degrees of freedom in the non-augmented problem.
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.
~FoldHandler()
Destructor, return the problem to its original state before the augmented system was added.
void solve_block_system()
Set to solve the block system.
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Get the residuals.
Vector< double > Y
Storage for the null vector.
void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tra...
void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable".
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...
int bifurcation_type() const
Indicate that we are tracking a fold bifurcation by returning 1.
FoldHandler(Problem *const &problem_pt, double *const &parameter_pt)
Constructor: initialise the fold handler, by setting initial guesses for Y, Phi and calculating count...
double * Parameter_pt
Storage for the pointer to the parameter.
unsigned ndof(GeneralisedElement *const &elem_pt)
Get the number of elemental degrees of freedom.
A Generalised Element class.
Definition: elements.h:73
A class that is used to assemble the augmented system that defines a Hopf bifurcation....
double * bifurcation_parameter_pt() const
Return a pointer to the bifurcation parameter in bifurcation tracking problems.
Vector< double > Phi
The real part of the null vector.
Vector< int > Count
A vector that is used to determine how many elements contribute to a particular equation....
double Omega
The critical frequency of the bifurcation.
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...
~HopfHandler()
Destructor, return the problem to its original state, before the augmented system was added.
void solve_full_system()
Solve non-block system.
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)
Get the residuals.
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.
int bifurcation_type() const
Indicate that we are tracking a Hopf bifurcation by returning 3.
unsigned ndof(GeneralisedElement *const &elem_pt)
Get the number of elemental degrees of freedom.
const double & omega() const
Return the frequency of the bifurcation.
Vector< double > C
A constant vector used to ensure that the null vector is not trivial.
void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tra...
Vector< double > Psi
The imaginary part of the null vector.
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Get the global equation number of the local unknown.
unsigned Solve_which_system
Integer flag to indicate which system should be assembled. There are three possibilities....
Problem * Problem_pt
Pointer to the problem.
void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable".
HopfHandler(Problem *const &problem_pt, double *const &parameter_pt)
Constructor.
unsigned Ndof
Store the number of degrees of freedom in the non-augmented problem.
void solve_complex_system()
Set to solve the complex system.
void solve_standard_system()
Set to solve the standard system.
double * Parameter_pt
Pointer to the parameter.
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
Base class for all linear solvers. This merely defines standard interfaces for linear solvers,...
Definition: linear_solver.h:68
An OomphLibError object which should be thrown when an run-time error is encountered....
A class that is used to assemble the residuals in parallel by overloading the get_all_vectors_and_mat...
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...
unsigned ndof(GeneralisedElement *const &elem_pt)
Use underlying assembly handler to return the number of degrees of freedom in the element elem_pt.
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" f...
AssemblyHandler * Assembly_handler_pt
The original assembly handler.
ParallelResidualsHandler(AssemblyHandler *const &assembly_handler_pt)
Constructor, set the original assembly handler.
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.
~ParallelResidualsHandler()
Empty virtual destructor.
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...
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" for elem_pt. Overloaded to return t...
unsigned ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt Pass through to the original assembly ...
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...
double * Parameter_pt
The value of the parameter.
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.
AssemblyHandler * Assembly_handler_pt
The original assembly handler.
ParameterDerivativeHandler(AssemblyHandler *const &assembly_handler_pt, double *const &parameter_pt)
Store the original assembly handler and parameter.
A class that is used to assemble the augmented system that defines a pitchfork (symmetry-breaking) bi...
LinearAlgebraDistribution * Dof_distribution_pt
Store the original dof distribution.
Problem * Problem_pt
Pointer to the problem.
void synchronise()
Function that is used to perform any synchronisation required during the solution.
void solve_full_system()
Solve non-block system.
DoubleVectorWithHaloEntries Y
Storage for the null vector.
double * bifurcation_parameter_pt() const
Return a pointer to the bifurcation parameter in bifurcation tracking problems.
AssemblyHandler * Assembly_handler_pt
Pointer to the underlying (original) assembly handler.
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Get the residuals.
double Sigma
A slack variable used to specify the amount of antisymmetry in the solution.
unsigned Solve_which_system
Integer flag to indicate which system should be assembled. There are three possibilities....
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...
~PitchForkHandler()
Destructor, return the problem to its original state, before the augmented system was added.
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.
DoubleVectorWithHaloEntries C
A constant vector used to ensure that the null vector is not trivial.
unsigned Ndof
Store the number of degrees of freedom in the non-augmented problem.
PitchForkHandler(Problem *const &problem_pt, AssemblyHandler *const &assembly_handler_pt, double *const &parameter_pt, const DoubleVector &symmetry_vector)
Constructor, initialise the systems.
void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tra...
DoubleVectorWithHaloEntries Psi
A constant vector that is specifies the symmetry being broken.
double * Parameter_pt
Storage for the pointer to the parameter.
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_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable".
int bifurcation_type() const
Indicate that we are tracking a pitchfork bifurcation by returning 2.
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Get the global equation number of the local unknown.
bool Distributed
Boolean to indicate whether the problem is distributed.
unsigned global_eqn_number(const unsigned &i)
Function that is used to return map the global equations using the simplistic numbering scheme into t...
Vector< unsigned > Global_eqn_number
A vector that is used to map the global equations to their actual location in a distributed problem.
void solve_augmented_block_system()
Set to solve the augmented block system.
void solve_block_system()
Set to solve the block system.
DoubleVectorWithHaloEntries Count
A vector that is used to determine how many elements contribute to a particular equation....
unsigned ndof(GeneralisedElement *const &elem_pt)
Get the number of elemental degrees of freedom.
LinearAlgebraDistribution * Augmented_dof_distribution_pt
The augmented distribution.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
//////////////////////////////////////////////////////////////////// ////////////////////////////////...