assembly_handler.cc
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 //OOOMPH-LIB includes
31 #include "assembly_handler.h"
32 #include "elements.h"
33 #include "problem.h"
34 #include "mesh.h"
35 
36 namespace oomph
37 {
38  ///////////////////////////////////////////////////////////////////////
39  // Non-inline functions for the AssemblyHandler class
40  //////////////////////////////////////////////////////////////////////
41 
42  //===================================================================
43  ///Get the number of elemental degrees of freedom. Direct call
44  ///to the function in the element.
45  //===================================================================
46  unsigned AssemblyHandler::ndof(GeneralisedElement* const &elem_pt)
47  {return elem_pt->ndof();}
48 
49  //===================================================================
50  ///Get the vector of dofs in the element elem_pt at time level t
51  ///Direct call to the function in the element.
52  //===================================================================
54  const unsigned &t, Vector<double> &dof)
55  {return elem_pt->dof_vector(t,dof);}
56 
57  //===================================================================
58  ///Get the vector of pointers to dofs in the element elem_pt
59  ///Direct call to the function in the element.
60  //===================================================================
62  Vector<double*> &dof_pt)
63  {return elem_pt->dof_pt_vector(dof_pt);}
64 
65  /// \short Return the t-th level of storage associated with the i-th
66  /// (local) dof stored in the problem
67  double &AssemblyHandler::local_problem_dof(Problem* const &problem_pt,
68  const unsigned &t,
69  const unsigned &i)
70  {return *(problem_pt->dof_pt(i)+t);}
71 
72 
73 
74  //==================================================================
75  ///Get the global equation number of the local unknown. Direct call
76  ///to the function in the element.
77  //==================================================================
78  unsigned long AssemblyHandler::eqn_number(GeneralisedElement* const &elem_pt,
79  const unsigned &ieqn_local)
80  {return elem_pt->eqn_number(ieqn_local);}
81 
82  //==================================================================
83  ///Get the residuals by calling the underlying element's residuals
84  ///directly.
85  //=================================================================
87  Vector<double> &residuals)
88  {elem_pt->get_residuals(residuals);}
89 
90  //=======================================================================
91  /// Calculate the elemental Jacobian matrix "d equation / d variable" by
92  /// calling the element's get_jacobian function.
93  //======================================================================
95  Vector<double> &residuals,
96  DenseMatrix<double> &jacobian)
97  {elem_pt->get_jacobian(residuals,jacobian);}
98 
99  //=======================================================================
100  /// Calculate all desired vectors and matrices that are required by
101  /// the element by calling the get_jacobian function
102  //=======================================================================
104  GeneralisedElement* const &elem_pt,
106  {
107  get_jacobian(elem_pt,vec[0],matrix[0]);
108  }
109 
110  //=======================================================================
111  /// \short Calculate the derivative of the residuals with respect to
112  /// a parameter, by calling the elemental function
113  //======================================================================
115  GeneralisedElement* const &elem_pt,
116  double* const &parameter_pt,
117  Vector<double> &dres_dparam)
118  {
119  elem_pt->get_dresiduals_dparameter(parameter_pt,dres_dparam);
120  }
121 
122 
123  //=====================================================================
124  /// \short Calculate the derivative of the residuals and jacobian
125  /// with respect to a parameter by calling the elemental function
126  //========================================================================
128  GeneralisedElement* const &elem_pt,
129  double* const &parameter_pt,
130  Vector<double> &dres_dparam,
131  DenseMatrix<double> &djac_dparam)
132  {
133  elem_pt->get_djacobian_dparameter(parameter_pt,dres_dparam,djac_dparam);
134  }
135 
136  /// \short Calculate the product of the Hessian (derivative of Jacobian with
137  /// respect to all variables) an eigenvector, Y, and
138  /// other specified vectors, C
139  /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}
140  /// At the moment the dof pointer is passed in to enable
141  /// easy calculation of finite difference default
143  GeneralisedElement* const &elem_pt,
144  Vector<double> const &Y,
145  DenseMatrix<double> const &C,
146  DenseMatrix<double> &product)
147  {
148  elem_pt->get_hessian_vector_products(Y,C,product);
149  }
150 
151 
152  //=======================================================================
153  /// Return the eigenfunction(s) associated with the bifurcation that
154  /// has been detected in bifurcation tracking problems. Default
155  /// Broken implementation
156  //=======================================================================
158  {
159  std::ostringstream error_stream;
160  error_stream
161  <<
162  "There is no bifurcation parameter associated with the current assembly handler.\n"
163  <<
164  "Eigenfunction are only calculated by the Fold, PitchFork and Hopf Handlers"
165  << "\n";
166 
167  throw OomphLibError(error_stream.str(),
168  OOMPH_CURRENT_FUNCTION,
169  OOMPH_EXCEPTION_LOCATION);
170  return 0;
171  }
172 
173 
174  //=======================================================================
175  /// Return the eigenfunction(s) associated with the bifurcation that
176  /// has been detected in bifurcation tracking problems. Default
177  /// Broken implementation
178  //=======================================================================
180  Vector<DoubleVector> &eigenfunction)
181  {
182  std::ostringstream error_stream;
183  error_stream
184  <<
185  "There is no eigenfunction associated with the current assembly handler.\n"
186  <<
187  "Eigenfunction are only calculated by the Fold, PitchFork and Hopf Handlers"
188  << "\n";
189 
190  throw OomphLibError(error_stream.str(),
191  OOMPH_CURRENT_FUNCTION,
192  OOMPH_EXCEPTION_LOCATION);
193  }
194 
195  ///==========================================================================
196  /// Compute the inner products of the given vector of pairs of
197  /// history values over the element. The values of the index in the pair
198  /// may be the same.
199  //===========================================================================
201  Vector<std::pair<unsigned,unsigned> >
202  const &history_index,
203  Vector<double> &inner_product)
204  {
205  elem_pt->get_inner_products(history_index,inner_product);
206  }
207 
208  //==========================================================================
209  /// Compute the vectors that when taken as a dot product with
210  /// other history values give the inner product over the element.
211  /// In other words if we call get_inner_product_vectors({0,1},output)
212  /// output[0] will be a vector such that dofs.output[0] gives the inner
213  /// product of current dofs with themselves.
214  //==========================================================================
216  GeneralisedElement* const &elem_pt, Vector<unsigned> const &history_index,
217  Vector<Vector<double> > &inner_product_vector)
218  {
219  elem_pt->get_inner_product_vectors(history_index,inner_product_vector);
220  }
221 
222 
223 
224  ///////////////////////////////////////////////////////////////////////
225  // Non-inline functions for the ExplicitTimeStepHandler class
226  //////////////////////////////////////////////////////////////////////
227 
228 
229  //===================================================================
230  ///Get the number of elemental degrees of freedom. Direct call
231  ///to the function in the element.
232  //===================================================================
234  {return elem_pt->ndof();}
235 
236  //==================================================================
237  ///Get the global equation number of the local unknown. Direct call
238  ///to the function in the element.
239  //==================================================================
240  unsigned long ExplicitTimeStepHandler::
241  eqn_number(GeneralisedElement* const &elem_pt,const unsigned &ieqn_local)
242  {return elem_pt->eqn_number(ieqn_local);}
243 
244  //==================================================================
245  ///Call the element's residuals
246  //=================================================================
248  GeneralisedElement* const &elem_pt,
249  Vector<double> &residuals)
250  {
251  elem_pt->get_residuals(residuals);
252  }
253 
254  //=======================================================================
255  ///Replace get jacobian with the call to get the mass matrix
256  //======================================================================
258  Vector<double> &residuals,
259  DenseMatrix<double> &jacobian)
260  {
261  elem_pt->get_mass_matrix(residuals,jacobian);
262  }
263 
264 
265  //=======================================================================
266  /// Calculate all desired vectors and matrices that are required by
267  /// the problem by calling those of the underlying element.
268  //=======================================================================
270  GeneralisedElement* const &elem_pt,
272  {
273 #ifdef PARANOID
274  //Check dimension
275  if(matrix.size() != 1)
276  {
277  throw OomphLibError(
278  "ExplicitTimeSteps should return one matrix",
279  OOMPH_CURRENT_FUNCTION,
280  OOMPH_EXCEPTION_LOCATION);
281  }
282 #endif
283  //Get the residuals and <mass matrix
284  elem_pt->get_mass_matrix(vec[0],matrix[0]);
285  }
286 
287 
288 
289  ///////////////////////////////////////////////////////////////////////
290  // Non-inline functions for the EigenProblemHandler class
291  //////////////////////////////////////////////////////////////////////
292 
293 
294  //===================================================================
295  ///Get the number of elemental degrees of freedom. Direct call
296  ///to the function in the element.
297  //===================================================================
299  {return elem_pt->ndof();}
300 
301  //==================================================================
302  ///Get the global equation number of the local unknown. Direct call
303  ///to the function in the element.
304  //==================================================================
305  unsigned long EigenProblemHandler::
306  eqn_number(GeneralisedElement* const &elem_pt,const unsigned &ieqn_local)
307  {return elem_pt->eqn_number(ieqn_local);}
308 
309  //==================================================================
310  ///Cannot call get_residuals for an eigenproblem, so throw an error
311  //=================================================================
313  Vector<double> &residuals)
314  {
315  throw OomphLibError(
316  "An eigenproblem does not have a get_residuals function",
317  OOMPH_CURRENT_FUNCTION,
318  OOMPH_EXCEPTION_LOCATION);
319  }
320 
321  //=======================================================================
322  ///Cannot call get_jacobian for an eigenproblem, so throw an error
323  //======================================================================
325  Vector<double> &residuals,
326  DenseMatrix<double> &jacobian)
327  {
328  throw OomphLibError(
329  "An eigenproblem does not have a get_jacobian function",
330  OOMPH_CURRENT_FUNCTION,
331  OOMPH_EXCEPTION_LOCATION);
332  }
333 
334 
335  //=======================================================================
336  /// Calculate all desired vectors and matrices that are required by
337  /// the problem by calling those of the underlying element.
338  //=======================================================================
340  GeneralisedElement* const &elem_pt,
342  {
343 #ifdef PARANOID
344  //Check dimension
345  if(matrix.size() != 2)
346  {
347  throw OomphLibError("EigenProblems should return two matrices",
348  OOMPH_CURRENT_FUNCTION,
349  OOMPH_EXCEPTION_LOCATION);
350  }
351 #endif
352  //Find the number of variables
353  unsigned n_var = elem_pt->ndof();
354  //Assign a dummy residuals vector
355  Vector<double> dummy(n_var);
356  //Get the jacobian and mass matrices
357  elem_pt->get_jacobian_and_mass_matrix(dummy,matrix[0],matrix[1]);
358 
359  //If we have a non-zero shift, then shift the A matrix
360  if(Sigma_real != 0.0)
361  {
362  //Set the shifted matrix
363  for(unsigned i=0;i<n_var;i++)
364  {
365  for(unsigned j=0;j<n_var;j++)
366  {
367  matrix[0](i,j) -= Sigma_real*matrix[1](i,j);
368  }
369  }
370  }
371  }
372 
373 
374  //======================================================================
375  /// Clean up the memory that may have been allocated by the solver
376  //=====================================================================
378  {
379  if(Alpha_pt!=0) {delete Alpha_pt;}
380  if(E_pt!=0) {delete E_pt;}
381  }
382 
383  //===================================================================
384  /// Use a block factorisation to solve the augmented system
385  /// associated with a PitchFork bifurcation.
386  //===================================================================
388  DoubleVector &result)
389  {
390  // if the result is setup then it should not be distributed
391 #ifdef PARANOID
392  if (result.built())
393  {
394  if (result.distributed())
395  {
396  throw OomphLibError("The result vector must not be distributed",
397  OOMPH_CURRENT_FUNCTION,
398  OOMPH_EXCEPTION_LOCATION);
399  }
400  }
401 #endif
402 
403  //Locally cache the pointer to the handler.
404  FoldHandler* handler_pt =
405  static_cast<FoldHandler*>(problem_pt->assembly_handler_pt());
406 
407  //Switch the handler to "block solver" mode
408  handler_pt->solve_augmented_block_system();
409 
410  //We need to find out the number of dofs in the problem
411  unsigned n_dof = problem_pt->ndof();
412 
413  // create the linear algebra distribution for this solver
414  // currently only global (non-distributed) distributions are allowed
415  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
416  this->build_distribution(dist);
417 
418  // if the result vector is not setup then rebuild with distribution = global
419  if (!result.built())
420  {
421  result.build(this->distribution_pt(),0.0);
422  }
423 
424  //Setup storage for temporary vectors
425  DoubleVector a(this->distribution_pt(),0.0), b(this->distribution_pt(),0.0);
426 
427  //Allocate storage for Alpha which can be used in the resolve
428  if(Alpha_pt!=0) {delete Alpha_pt;}
429  Alpha_pt = new DoubleVector(this->distribution_pt(),0.0);
430 
431  //We are going to do resolves using the underlying linear solver
433 
434  //Solve the first system Aa = R
435  Linear_solver_pt->solve(problem_pt,a);
436 
437  //The vector in the top right-hand block is the jacobian multiplied
438  //by the null vector
439 
440  //Get the present null vector from the handler
441  DoubleVector y(this->distribution_pt(),0.0);
442  for(unsigned n=0;n<(n_dof-1);++n) {y[n] = handler_pt->Y[n];}
443  //For simplicity later, add a zero at the end
444  y[n_dof-1] = 0.0;
445 
446  //Loop over the elements and assemble the product
447  //Local storage for the product terms
448  DoubleVector Jy(this->distribution_pt(),0.0);
449 
450  //Calculate the product of the jacobian matrices, etc
451  unsigned long n_element = problem_pt->mesh_pt()->nelement();
452  for(unsigned long e = 0;e<n_element;e++)
453  {
454  GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
455  //Loop over the ndofs in each element
456  unsigned n_var = elem_pt->ndof();
457  //Get the jacobian matrix
458  DenseMatrix<double> jac(n_var);
459  // storage for the residual
460  Vector<double> res(n_var);
461  //Get unperturbed raw jacobian
462  elem_pt->get_jacobian(res,jac);
463 
464  //Multiply the dofs
465  for(unsigned n=0;n<n_var;n++)
466  {
467  unsigned eqn_number = elem_pt->eqn_number(n);
468  for(unsigned m=0;m<n_var;m++)
469  {
470  unsigned unknown = elem_pt->eqn_number(m);
471  Jy[eqn_number] += jac(n,m)*y[unknown];
472  }
473  }
474  }
475  //The final entry of the vector will be zero
476 
477  //Now resolve to find alpha
479 
480  //The vector that multiplie the product matrix is actually y - alpha
481  DoubleVector y_minus_alpha(this->distribution_pt(),0.0);
482  for(unsigned n=0;n<n_dof;n++)
483  {
484  y_minus_alpha[n] = y[n] - (*Alpha_pt)[n];
485  }
486 
487  //We can now construct our multipliers
488  //Prepare to scale
489  double dof_length=0.0, a_length=0.0, alpha_length=0.0;
490  for(unsigned n=0;n<n_dof;n++)
491  {
492  if(std::fabs(problem_pt->dof(n)) > dof_length)
493  {dof_length = std::fabs(problem_pt->dof(n));}
494  if(std::fabs(a[n]) > a_length) {a_length = std::fabs(a[n]);}
495  if(std::fabs(y_minus_alpha[n]) > alpha_length)
496  {alpha_length = std::fabs(y_minus_alpha[n]);}
497  }
498 
499  double a_mult = dof_length/a_length;
500  double alpha_mult = dof_length/alpha_length;
501  const double FD_step = 1.0e-8;
502  a_mult += FD_step; alpha_mult += FD_step;
503  a_mult *= FD_step; alpha_mult *= FD_step;
504 
505  //Local storage for the product terms
506  DoubleVector Jprod_a(this->distribution_pt(),0.0),
507  Jprod_alpha(this->distribution_pt(),0.0);
508 
509  //Calculate the product of the jacobian matrices, etc
510  for(unsigned long e = 0;e<n_element;e++)
511  {
512  GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
513  //Loop over the ndofs in each element
514  unsigned n_var = handler_pt->ndof(elem_pt);
515  //Get the jacobian matrices
516  DenseMatrix<double> jac(n_var), jac_a(n_var), jac_alpha(n_var);
517  // elemental residual storage
518  Vector<double> res_elemental(n_var);
519  //Get unperturbed jacobian
520  handler_pt->get_jacobian(elem_pt,res_elemental,jac);
521 
522  //Backup the dofs
523  Vector<double> dof_bac(n_var);
524  //Perturb the dofs
525  for(unsigned n=0;n<n_var;n++)
526  {
527  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
528  dof_bac[n] = problem_pt->dof(eqn_number);
529  //Pertub by vector a
530  problem_pt->dof(eqn_number) += a_mult*a[eqn_number];
531  }
532 
534 
535  //Now get the new jacobian
536  handler_pt->get_jacobian(elem_pt,res_elemental,jac_a);
537 
538  //Perturb the dofs
539  for(unsigned n=0;n<n_var;n++)
540  {
541  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
542  problem_pt->dof(eqn_number) = dof_bac[n];
543  //Pertub by vector a
544  problem_pt->dof(eqn_number) += alpha_mult*y_minus_alpha[eqn_number];
545  }
546 
548 
549  //Now get the new jacobian
550  handler_pt->get_jacobian(elem_pt,res_elemental,jac_alpha);
551 
552  //Reset the dofs
553  for(unsigned n=0;n<n_var;n++)
554  {
555  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
556  problem_pt->dof(eqn_number) = dof_bac[n];
557  }
558 
560 
561  //OK, now work out the products
562  //Note the (n_var-1), we are only interested in the non-augmented
563  //jacobian
564  for(unsigned n=0;n<(n_var-1);n++)
565  {
566  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
567  double prod_a=0.0, prod_alpha=0.0;
568  for(unsigned m=0;m<(n_var-1);m++)
569  {
570  unsigned unknown = handler_pt->eqn_number(elem_pt,m);
571  prod_a += (jac_a(n,m) - jac(n,m))*y[unknown];
572  prod_alpha += (jac_alpha(n,m) - jac(n,m))*y[unknown];
573  }
574  Jprod_a[eqn_number] += prod_a/a_mult;
575  Jprod_alpha[eqn_number] += prod_alpha/alpha_mult;
576  }
577  }
578 
579  Jprod_alpha[n_dof-1] = 0.0;
580  Jprod_a[n_dof-1] = 0.0;
581 
582  //OK, now we can formulate the next vectors
583  //The assumption here is that the result has been set to the
584  //residuals.
585  for(unsigned n=0;n<n_dof-1;n++)
586  {
587  b[n] = result[n_dof+n] - Jprod_a[n];
588  }
589  //The final residual is the entry corresponding to the original parameter
590  b[n_dof-1] = result[n_dof-1];
591 
592  //Allocate storage for E which can be used in the resolve
593  if(E_pt!=0) {delete E_pt;}
594  E_pt = new DoubleVector(this->distribution_pt(),0.0);
595  DoubleVector f(this->distribution_pt(),0.0);
597  Linear_solver_pt->resolve(Jprod_alpha,*E_pt);
598 
599  //Calculate the final entry in the vector e
600  const double e_final = (*E_pt)[n_dof-1];
601  //Calculate the final entry in the vector d
602  const double d_final = f[n_dof-1]/e_final;
603  //Assemble the final corrections
604  for(unsigned n=0;n<n_dof-1;n++)
605  {
606  result[n] = a[n] - (*Alpha_pt)[n]*d_final + d_final*y[n];
607  result[n_dof+n] = f[n] - (*E_pt)[n]*d_final;
608  }
609  //The result corresponding to the parameter
610  result[n_dof-1] = a[n_dof-1] - (*Alpha_pt)[n_dof-1]*d_final;
611 
612  //The sign of the jacobian is the sign of the final entry in e
613  problem_pt->sign_of_jacobian() =
614  static_cast<int>(std::fabs(e_final)/e_final);
615 
616  //Switch things to our block solver
617  handler_pt->solve_full_system();
618 
619  //If we are not storing things, clear the results
620  if(!Enable_resolve)
621  {
622  //We no longer need to store the matrix
624  delete Alpha_pt; Alpha_pt=0;
625  delete E_pt; E_pt=0;
626  }
627  //Otherwise, also store the problem pointer
628  else
629  {
630  Problem_pt = problem_pt;
631  }
632  }
633 
634 
635  //======================================================================
636  //Hack the re-solve to use the block factorisation
637  //======================================================================
639  DoubleVector &result)
640  {
641  //Check that the factors have been stored
642  if(Alpha_pt==0)
643  {
644  throw OomphLibError("The required vectors have not been stored",
645  OOMPH_CURRENT_FUNCTION,
646  OOMPH_EXCEPTION_LOCATION);
647  }
648 
649  //Get the pointer to the problem
650  Problem* const problem_pt = Problem_pt;
651 
652  FoldHandler* handler_pt =
653  static_cast<FoldHandler*>(problem_pt->assembly_handler_pt());
654 
655  //Switch things to our block solver
656  handler_pt->solve_augmented_block_system();
657 
658  //We need to find out the number of dofs in the problem
659  unsigned n_dof = problem_pt->ndof();
660 
661 #ifdef PARANOID
662  // if the result is setup then it should not be distributed
663  if (result.built())
664  {
665  if (result.distributed())
666  {
667  throw OomphLibError("The result vector must not be distributed",
668  OOMPH_CURRENT_FUNCTION,
669  OOMPH_EXCEPTION_LOCATION);
670  }
671  }
672  // the rhs must be setup
673  if (!rhs.built())
674  {
675  throw OomphLibError("The rhs vector must be setup",
676  OOMPH_CURRENT_FUNCTION,
677  OOMPH_EXCEPTION_LOCATION);
678  }
679 #endif
680 
681  // if the result vector is not setup then rebuild with distribution = global
682  if (!result.built())
683  {
684  result.build(rhs.distribution_pt(),0.0);
685  }
686 
687  //Setup storage
688  DoubleVector a(this->distribution_pt(),0.0), b(this->distribution_pt(),0.0);
689 
690  //Set the values of the a vector
691  for(unsigned n=0;n<(n_dof-1);n++) {a[n] = rhs[n];}
692  //The entry associated with the additional parameter is zero
693  a[n_dof-1] = 0.0;
694 
696 
697  // Copy rhs vector into local storage so it doesn't get overwritten
698  // if the linear solver decides to initialise the solution vector, say,
699  // which it's quite entitled to do!
700  DoubleVector input_a(a);
701 
702  Linear_solver_pt->resolve(input_a,a);
703 
704  //We can now construct our multipliers
705  //Prepare to scale
706  double dof_length=0.0, a_length=0.0;
707  for(unsigned n=0;n<n_dof;n++)
708  {
709  if(std::fabs(problem_pt->dof(n)) > dof_length)
710  {dof_length = std::fabs(problem_pt->dof(n));}
711 
712  if(std::fabs(a[n]) > a_length) {a_length = std::fabs(a[n]);}
713  }
714  double a_mult = dof_length/a_length;
715  const double FD_step = 1.0e-8;
716  a_mult += FD_step;
717  a_mult *= FD_step;
718 
719  DoubleVector Jprod_a(this->distribution_pt(),0.0);
720 
721  unsigned long n_element = problem_pt->mesh_pt()->nelement();
722  for(unsigned long e = 0;e<n_element;e++)
723  {
724  GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
725  //Loop over the ndofs in each element
726  unsigned n_var = handler_pt->ndof(elem_pt);
727  //Get some jacobian matrices
728  DenseMatrix<double> jac(n_var), jac_a(n_var);
729  // the elemental residual
730  Vector<double> res_elemental(n_var);
731  //Get unperturbed jacobian
732  handler_pt->get_jacobian(elem_pt,res_elemental,jac);
733 
734  //Backup the dofs
735  Vector<double> dof_bac(n_var);
736  //Perturb the dofs
737  for(unsigned n=0;n<n_var;n++)
738  {
739  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
740  dof_bac[n] = problem_pt->dof(eqn_number);
741  //Pertub by vector a
742  problem_pt->dof(eqn_number) += a_mult*a[eqn_number];
743  }
744 
746 
747  //Now get the new jacobian
748  handler_pt->get_jacobian(elem_pt,res_elemental,jac_a);
749 
750  //Reset the dofs
751  for(unsigned n=0;n<n_var;n++)
752  {
753  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
754  problem_pt->dof(eqn_number) = dof_bac[n];
755  }
756 
758 
759  //OK, now work out the products
760  for(unsigned n=0;n<(n_var-1);n++)
761  {
762  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
763  double prod_a=0.0;
764  for(unsigned m=0;m<(n_var-1);m++)
765  {
766  unsigned unknown = handler_pt->eqn_number(elem_pt,m);
767  prod_a += (jac_a(n,m) - jac(n,m))*handler_pt->Y[unknown];
768  }
769  Jprod_a[eqn_number] += prod_a/a_mult;
770  }
771  }
772 
773  Jprod_a[n_dof-1] = 0.0;
774 
775  //OK, now we can formulate the next vectors
776  for(unsigned n=0;n<n_dof-1;n++)
777  {
778  b[n] = rhs[n_dof+n] - Jprod_a[n];
779  }
780  //The residuals for the final entry should be those associated
781  //with the parameter
782  b[n_dof-1] = rhs[n_dof-1];
783 
784  DoubleVector f(this->distribution_pt(),0.0);
785 
787 
788  //Calculate the final entry in the vector d
789  const double d_final = f[n_dof-1]/(*E_pt)[n_dof-1];
790  //Assemble the final corrections
791  for(unsigned n=0;n<n_dof-1;n++)
792  {
793  result[n] = a[n] - (*Alpha_pt)[n]*d_final + d_final*handler_pt->Y[n];
794  result[n_dof+n] = f[n] - (*E_pt)[n]*d_final;
795  }
796  //The result corresponding to the paramater
797  result[n_dof-1] = a[n_dof-1] - (*Alpha_pt)[n_dof-1]*d_final;
798 
800 
801  //Switch things to our block solver
802  handler_pt->solve_full_system();
803  }
804 
805 
806  ///////////////////////////////////////////////////////////////////////
807  // Non-inline functions for the FoldHandler class
808  //////////////////////////////////////////////////////////////////////
809 
810 
811  //========================================================================
812  /// Constructor: Initialise the fold handler,
813  /// by setting initial guesses for Y, Phi and calculating Count.
814  /// If the system changes, a new handler must be constructed.
815  //========================================================================
816  FoldHandler::FoldHandler(Problem* const &problem_pt,
817  double* const &parameter_pt) :
818  Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
819  {
820  //Set the problem pointer
821  Problem_pt = problem_pt;
822  //Set the number of degrees of freedom
823  Ndof = problem_pt->ndof();
824 
825  // create the linear algebra distribution for this solver
826  // currently only global (non-distributed) distributions are allowed
827  LinearAlgebraDistribution* dist_pt = new
828  LinearAlgebraDistribution(problem_pt->communicator_pt(),Ndof,false);
829 
830  //Resize the vectors of additional dofs and constants
831  Phi.resize(Ndof);
832  Y.resize(Ndof);
833  Count.resize(Ndof,0);
834 
835  //Loop over all the elements in the problem
836  unsigned n_element = problem_pt->mesh_pt()->nelement();
837  for(unsigned e=0;e<n_element;e++)
838  {
839  GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
840  //Loop over the local freedoms in the element
841  unsigned n_var = elem_pt->ndof();
842  for(unsigned n=0;n<n_var;n++)
843  {
844  //Increase the associated global equation number counter
845  ++Count[elem_pt->eqn_number(n)];
846  }
847  }
848 
849  //Calculate the value Phi by
850  //solving the system JPhi = dF/dlambda
851 
852  //Locally cache the linear solver
853  LinearSolver* const linear_solver_pt = problem_pt->linear_solver_pt();
854 
855  //Save the status before entry to this routine
856  bool enable_resolve = linear_solver_pt->is_resolve_enabled();
857 
858  //We need to do a resolve
859  linear_solver_pt->enable_resolve();
860 
861  //Storage for the solution
862  DoubleVector x(dist_pt,0.0);
863 
864  //Solve the standard problem, we only want to make sure that
865  //we factorise the matrix, if it has not been factorised. We shall
866  //ignore the return value of x.
867  linear_solver_pt->solve(problem_pt,x);
868 
869  //Get the vector dresiduals/dparameter
870  problem_pt->get_derivative_wrt_global_parameter(parameter_pt,x);
871 
872  // Copy rhs vector into local storage so it doesn't get overwritten
873  // if the linear solver decides to initialise the solution vector, say,
874  // which it's quite entitled to do!
875  DoubleVector input_x(x);
876 
877  //Now resolve the system with the new RHS and overwrite the solution
878  linear_solver_pt->resolve(input_x,x);
879 
880  //Restore the storage status of the linear solver
881  if (enable_resolve)
882  {
883  linear_solver_pt->enable_resolve();
884  }
885  else
886  {
887  linear_solver_pt->disable_resolve();
888  }
889 
890  //Add the global parameter as an unknown to the problem
891  problem_pt->Dof_pt.push_back(parameter_pt);
892 
893 
894  //Normalise the initial guesses for phi
895  double length = 0.0;
896  for(unsigned n=0;n<Ndof;n++) {length += x[n]*x[n];}
897  length = sqrt(length);
898 
899  //Now add the null space components to the problem unknowns
900  //and initialise them and Phi to the same normalised values
901  for(unsigned n=0;n<Ndof;n++)
902  {
903  problem_pt->Dof_pt.push_back(&Y[n]);
904  Y[n] = Phi[n] = -x[n]/length;
905  }
906 
907  // delete the dist_pt
908  problem_pt->Dof_distribution_pt->build(problem_pt->communicator_pt(),
909  Ndof*2+1,true);
910  //Remove all previous sparse storage used during Jacobian assembly
912 
913  delete dist_pt;
914  }
915 
916 
917  ///Constructor in which the eigenvector is passed as an initial
918  ///guess
919  FoldHandler::FoldHandler(Problem* const &problem_pt,
920  double* const &parameter_pt,
921  const DoubleVector &eigenvector) :
922  Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
923  {
924  //Set the problem pointer
925  Problem_pt = problem_pt;
926  //Set the number of degrees of freedom
927  Ndof = problem_pt->ndof();
928 
929  // create the linear algebra distribution for this solver
930  // currently only global (non-distributed) distributions are allowed
931  LinearAlgebraDistribution* dist_pt = new
933  Ndof,false);
934 
935  //Resize the vectors of additional dofs and constants
936  Phi.resize(Ndof);
937  Y.resize(Ndof);
938  Count.resize(Ndof,0);
939 
940  //Loop over all the elements in the problem
941  unsigned n_element = problem_pt->mesh_pt()->nelement();
942  for(unsigned e=0;e<n_element;e++)
943  {
944  GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
945  //Loop over the local freedoms in the element
946  unsigned n_var = elem_pt->ndof();
947  for(unsigned n=0;n<n_var;n++)
948  {
949  //Increase the associated global equation number counter
950  ++Count[elem_pt->eqn_number(n)];
951  }
952  }
953 
954  //Add the global parameter as an unknown to the problem
955  problem_pt->Dof_pt.push_back(parameter_pt);
956 
957 
958  //Normalise the initial guesses for the eigenvecor
959  double length = 0.0;
960  for(unsigned n=0;n<Ndof;n++)
961  {length += eigenvector[n]*eigenvector[n];}
962  length = sqrt(length);
963 
964  //Now add the null space components to the problem unknowns
965  //and initialise them and Phi to the same normalised values
966  for(unsigned n=0;n<Ndof;n++)
967  {
968  problem_pt->Dof_pt.push_back(&Y[n]);
969  Y[n] = Phi[n] = eigenvector[n]/length;
970  }
971 
972  // delete the dist_pt
973  problem_pt->Dof_distribution_pt->build(problem_pt->communicator_pt(),
974  Ndof*2+1,true);
975  //Remove all previous sparse storage used during Jacobian assembly
977 
978  delete dist_pt;
979  }
980 
981 
982  ///Constructor in which the eigenvector and normalisation
983  ///vector are passed as an initial guess
984  FoldHandler::FoldHandler(Problem* const &problem_pt,
985  double* const &parameter_pt,
986  const DoubleVector &eigenvector,
987  const DoubleVector &normalisation) :
988  Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
989  {
990  //Set the problem pointer
991  Problem_pt = problem_pt;
992  //Set the number of degrees of freedom
993  Ndof = problem_pt->ndof();
994 
995  // create the linear algebra distribution for this solver
996  // currently only global (non-distributed) distributions are allowed
997  LinearAlgebraDistribution* dist_pt = new
999  Ndof,false);
1000 
1001  //Resize the vectors of additional dofs and constants
1002  Phi.resize(Ndof);
1003  Y.resize(Ndof);
1004  Count.resize(Ndof,0);
1005 
1006  //Loop over all the elements in the problem
1007  unsigned n_element = problem_pt->mesh_pt()->nelement();
1008  for(unsigned e=0;e<n_element;e++)
1009  {
1010  GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
1011  //Loop over the local freedoms in the element
1012  unsigned n_var = elem_pt->ndof();
1013  for(unsigned n=0;n<n_var;n++)
1014  {
1015  //Increase the associated global equation number counter
1016  ++Count[elem_pt->eqn_number(n)];
1017  }
1018  }
1019 
1020  //Add the global parameter as an unknown to the problem
1021  problem_pt->Dof_pt.push_back(parameter_pt);
1022 
1023 
1024  //Normalise the initial guesses for the eigenvecor
1025  double length = 0.0;
1026  for(unsigned n=0;n<Ndof;n++)
1027  {length += eigenvector[n]*normalisation[n];}
1028  length = sqrt(length);
1029 
1030  //Now add the null space components to the problem unknowns
1031  //and initialise them and Phi to the same normalised values
1032  for(unsigned n=0;n<Ndof;n++)
1033  {
1034  problem_pt->Dof_pt.push_back(&Y[n]);
1035  Y[n] = eigenvector[n]/length;
1036  Phi[n] = normalisation[n];
1037  }
1038 
1039  // delete the dist_pt
1040  problem_pt->Dof_distribution_pt->build(problem_pt->communicator_pt(),
1041  Ndof*2+1,true);
1042  //Remove all previous sparse storage used during Jacobian assembly
1044 
1045  delete dist_pt;
1046  }
1047 
1048 
1049 
1050  //=================================================================
1051  /// The number of unknowns is 2n+1
1052  //================================================================
1053  unsigned FoldHandler::ndof(GeneralisedElement* const &elem_pt)
1054  {
1055  unsigned raw_ndof = elem_pt->ndof();
1056  //Return different values depending on the type of block decomposition
1057  switch(Solve_which_system)
1058  {
1059  case Full_augmented:
1060  return (2*raw_ndof + 1);
1061  break;
1062 
1063  case Block_augmented_J:
1064  return (raw_ndof+1);
1065  break;
1066 
1067  case Block_J:
1068  return raw_ndof;
1069  break;
1070 
1071  default:
1072  std::ostringstream error_stream;
1073  error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
1074  << " not " << Solve_which_system << "\n";
1075  throw OomphLibError(error_stream.str(),
1076  OOMPH_CURRENT_FUNCTION,
1077  OOMPH_EXCEPTION_LOCATION);
1078  }
1079  }
1080 
1081  //=====================================================================
1082  /// Return the global equation number associated with local equation
1083  /// number ieqn_local. We choose to number the unknowns according
1084  /// to the augmented system.
1085  //=======================================================================
1086  unsigned long FoldHandler::eqn_number(GeneralisedElement* const &elem_pt,
1087  const unsigned &ieqn_local)
1088  {
1089  //Find the number of non-augmented dofs in the element
1090  unsigned raw_ndof = elem_pt->ndof();
1091  //Storage for the global eqn number
1092  unsigned long global_eqn=0;
1093  //If we are a "standard" unknown, just return the standard equation number
1094  if(ieqn_local < raw_ndof)
1095  {
1096  global_eqn = elem_pt->eqn_number(ieqn_local);
1097  }
1098  //Otherwise if we are at an unknown corresponding to the bifurcation
1099  //parameter return the global equation number of the parameter
1100  else if(ieqn_local==raw_ndof)
1101  {
1102  global_eqn = Ndof;
1103  }
1104  //Otherwise we are in the unknown corresponding to a null vector
1105  //return the global unknown Ndof + 1 + global unknown of "standard" unknown.
1106  else
1107  {
1108  global_eqn = Ndof + 1 + elem_pt->eqn_number(ieqn_local - 1 - raw_ndof);
1109  }
1110 
1111  //Return the global equation number
1112  return global_eqn;
1113  }
1114 
1115  //====================================================================
1116  ///Formulate the augmented system
1117  //====================================================================
1119  Vector<double> &residuals)
1120  {
1121  //Need to get raw residuals and jacobian
1122  unsigned raw_ndof = elem_pt->ndof();
1123 
1124  //Find out which system we are solving
1125  switch(Solve_which_system)
1126  {
1127  //If we are solving the standard system
1128  case Block_J:
1129  {
1130  //Get the basic residuals
1131  elem_pt->get_residuals(residuals);
1132  }
1133  break;
1134 
1135  //If we are solving the augmented-by-one system
1136  case Block_augmented_J:
1137  {
1138  //Get the basic residuals
1139  elem_pt->get_residuals(residuals);
1140 
1141  //Zero the final residual
1142  residuals[raw_ndof] = 0.0;
1143  }
1144  break;
1145 
1146  //If we are solving the full augmented system
1147  case Full_augmented:
1148  {
1149  DenseMatrix<double> jacobian(raw_ndof);
1150  //Get the basic residuals and jacobian initially
1151  elem_pt->get_jacobian(residuals,jacobian);
1152 
1153  //The normalisation equation must be initialised to -1.0/number of elements
1154  residuals[raw_ndof] = -1.0/Problem_pt->mesh_pt()->nelement();
1155 
1156  //Now assemble the equations Jy = 0
1157  for(unsigned i=0;i<raw_ndof;i++)
1158  {
1159  residuals[raw_ndof+1+i] = 0.0;
1160  for(unsigned j=0;j<raw_ndof;j++)
1161  {
1162  residuals[raw_ndof+1+i] += jacobian(i,j)*Y[elem_pt->eqn_number(j)];
1163  }
1164  //Add the contribution to phi.y=1
1165  //Need to divide by the number of elements that contribute to this
1166  //unknown so that we do get phi.y exactly.
1167  unsigned global_eqn = elem_pt->eqn_number(i);
1168  residuals[raw_ndof] +=
1169  (Phi[global_eqn]*Y[global_eqn])/Count[global_eqn];
1170  }
1171  }
1172  break;
1173 
1174  default:
1175  std::ostringstream error_stream;
1176  error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
1177  << " not " << Solve_which_system << "\n";
1178  throw OomphLibError(error_stream.str(),
1179  OOMPH_CURRENT_FUNCTION,
1180  OOMPH_EXCEPTION_LOCATION);
1181  }
1182  }
1183 
1184 
1185  //=================================================================
1186  /// Calculate the elemental Jacobian matrix "d equation
1187  /// / d variable" for the augmented system
1188  //=================================================================
1190  Vector<double> &residuals,
1191  DenseMatrix<double> &jacobian)
1192  {
1193 
1194  //Find the number of augmented dofs
1195  unsigned augmented_ndof = ndof(elem_pt);
1196  //Find the non-augmented dofs
1197  unsigned raw_ndof = elem_pt->ndof();
1198 
1199  //Which system are we solving
1200  switch(Solve_which_system)
1201  {
1202  //If we are solving the original system
1203  case Block_J:
1204  {
1205  //Just get the raw jacobian and residuals
1206  elem_pt->get_jacobian(residuals,jacobian);
1207  }
1208  break;
1209 
1210  //If we are solving the augmented-by-one system
1211  case Block_augmented_J:
1212  {
1213  //Get the full residuals, we need them
1214  get_residuals(elem_pt,residuals);
1215 
1216  //Need to get the raw jacobian (and raw residuals)
1217  Vector<double> newres(augmented_ndof);
1218  elem_pt->get_jacobian(newres,jacobian);
1219 
1220  //Now do finite differencing stuff
1221  const double FD_step = 1.0e-8;
1222  //Fill in the first lot of finite differences
1223  {
1224  //increase the global parameter
1225  double *unknown_pt = Problem_pt->Dof_pt[Ndof];
1226  double init = *unknown_pt;
1227  *unknown_pt += FD_step;
1228 
1229  //Now do any possible updates
1230  Problem_pt
1232 
1233  //Get the new (modified) residuals
1234  get_residuals(elem_pt,newres);
1235 
1236  //The final column is given by the difference
1237  //between the residuals
1238  for(unsigned n=0;n<raw_ndof;n++)
1239  {
1240  jacobian(n,augmented_ndof-1) = (newres[n] - residuals[n])/FD_step;
1241  }
1242  //Reset the global parameter
1243  *unknown_pt = init;
1244 
1245  //Now do any possible updates
1246  Problem_pt
1248  }
1249 
1250  //Fill in the bottom row
1251  for(unsigned n=0;n<raw_ndof;n++)
1252  {
1253  unsigned local_eqn = elem_pt->eqn_number(n);
1254  jacobian(augmented_ndof-1,n)
1255  = Phi[local_eqn]/Count[local_eqn];
1256  }
1257  }
1258  break;
1259 
1260  //Otherwise solving the full system
1261  case Full_augmented:
1262  {
1263  //Get the residuals for the augmented system
1264  get_residuals(elem_pt,residuals);
1265 
1266  //Need to get the raw residuals and jacobian
1267  Vector<double> newres(raw_ndof);
1268  DenseMatrix<double> newjac(raw_ndof);
1269  elem_pt->get_jacobian(newres,jacobian);
1270 
1271  //Fill in the jacobian on the diagonal sub-block of
1272  //the null-space equations
1273  for(unsigned n=0;n<raw_ndof;n++)
1274  {
1275  for(unsigned m=0;m<raw_ndof;m++)
1276  {
1277  jacobian(raw_ndof+1+n,raw_ndof+1+m) = jacobian(n,m);
1278  }
1279  }
1280 
1281  //Now finite difference wrt the global unknown
1282  const double FD_step = 1.0e-8;
1283  //Fill in the first lot of finite differences
1284  {
1285  //increase the global parameter
1286  double *unknown_pt = Problem_pt->Dof_pt[Ndof];
1287  double init = *unknown_pt;
1288  *unknown_pt += FD_step;
1289  //Need to update the function
1291 
1292  //Get the new raw residuals and jacobian
1293  elem_pt->get_jacobian(newres,newjac);
1294 
1295  //The end of the first row is given by the difference
1296  //between the residuals
1297  for(unsigned n=0;n<raw_ndof;n++)
1298  {
1299  jacobian(n,raw_ndof) = (newres[n] - residuals[n])/FD_step;
1300  //The end of the second row is given by the difference multiplied by
1301  //the product
1302  for(unsigned l=0;l<raw_ndof;l++)
1303  {
1304  jacobian(raw_ndof+1+n,raw_ndof)
1305  += (newjac(n,l) - jacobian(n,l))*Y[elem_pt->eqn_number(l)]/
1306  FD_step;
1307  }
1308  }
1309  //Reset the global parameter
1310  *unknown_pt = init;
1311 
1312  //Need to update the function
1314  }
1315 
1316  //Now fill in the first column of the second rows
1317  {
1318  for(unsigned n=0;n<raw_ndof;n++)
1319  {
1320  unsigned long global_eqn = eqn_number(elem_pt,n);
1321  //Increase the first lot
1322  double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
1323  double init = *unknown_pt;
1324  *unknown_pt += FD_step;
1326 
1327  //Get the new jacobian
1328  elem_pt->get_jacobian(newres,newjac);
1329 
1330  //Work out the differences
1331  for(unsigned k=0;k<raw_ndof;k++)
1332  {
1333  //jacobian(raw_ndof+k,n) = 0.0;
1334  for(unsigned l=0;l<raw_ndof;l++)
1335  {
1336  jacobian(raw_ndof+1+k,n)
1337  += (newjac(k,l) - jacobian(k,l))*Y[elem_pt->eqn_number(l)]/
1338  FD_step;
1339  }
1340  }
1341  *unknown_pt = init;
1343  }
1344  }
1345 
1346  //Fill in the row corresponding to the parameter
1347  for(unsigned n=0;n<raw_ndof;n++)
1348  {
1349  unsigned global_eqn = elem_pt->eqn_number(n);
1350  jacobian(raw_ndof,raw_ndof+1+n)
1351  = Phi[global_eqn]/Count[global_eqn];
1352  }
1353 
1354  // //Loop over the dofs
1355 // for(unsigned n=0;n<augmented_ndof;n++)
1356 // {
1357 // unsigned long global_eqn = eqn_number(elem_pt,n);
1358 // double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
1359 // double init = *unknown_pt;
1360 // *unknown_pt += FD_step;
1361 
1362 // //Get the new residuals
1363 // get_residuals(elem_pt,newres);
1364 
1365 // for(unsigned m=0;m<augmented_ndof;m++)
1366 // {
1367 // jacobian(m,n) = (newres[m] - residuals[m])/FD_step;
1368 // }
1369 // //Reset the unknown
1370 // *unknown_pt = init;
1371 // }
1372  }
1373  break;
1374 
1375  default:
1376  std::ostringstream error_stream;
1377  error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
1378  << " not " << Solve_which_system << "\n";
1379  throw OomphLibError(error_stream.str(),
1380  OOMPH_CURRENT_FUNCTION,
1381  OOMPH_EXCEPTION_LOCATION);
1382 
1383  }
1384  }
1385 
1386 
1387 
1388  //====================================================================
1389  ///Formulate the derivatives of the augmented system with respect
1390  ///to a parameter
1391  //====================================================================
1393  GeneralisedElement* const &elem_pt,
1394  double* const &parameter_pt,
1395  Vector<double> &dres_dparam)
1396  {
1397  //Need to get raw residuals and jacobian
1398  unsigned raw_ndof = elem_pt->ndof();
1399 
1400  //Find out which system we are solving
1401  switch(Solve_which_system)
1402  {
1403  //If we are solving the standard system
1404  case Block_J:
1405  {
1406  //Get the basic residual derivatives
1407  elem_pt->get_dresiduals_dparameter(parameter_pt,dres_dparam);
1408  }
1409  break;
1410 
1411  //If we are solving the augmented-by-one system
1412  case Block_augmented_J:
1413  {
1414  //Get the basic residual derivatives
1415  elem_pt->get_dresiduals_dparameter(parameter_pt,dres_dparam);
1416 
1417  //Zero the final derivative
1418  dres_dparam[raw_ndof] = 0.0;
1419  }
1420  break;
1421 
1422  //If we are solving the full augmented system
1423  case Full_augmented:
1424  {
1425  DenseMatrix<double> djac_dparam(raw_ndof);
1426  //Get the basic residuals and jacobian derivatives initially
1427  elem_pt->get_djacobian_dparameter(parameter_pt,dres_dparam,djac_dparam);
1428 
1429  //The normalisation equation does not depend on the parameter
1430  dres_dparam[raw_ndof] = 0.0;
1431 
1432  //Now assemble the equations dJy/dparameter = 0
1433  for(unsigned i=0;i<raw_ndof;i++)
1434  {
1435  dres_dparam[raw_ndof+1+i] = 0.0;
1436  for(unsigned j=0;j<raw_ndof;j++)
1437  {
1438  dres_dparam[raw_ndof+1+i] += djac_dparam(i,j)*
1439  Y[elem_pt->eqn_number(j)];
1440  }
1441  }
1442  }
1443  break;
1444 
1445  default:
1446  std::ostringstream error_stream;
1447  error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
1448  << " not " << Solve_which_system << "\n";
1449  throw OomphLibError(error_stream.str(),
1450  OOMPH_CURRENT_FUNCTION,
1451  OOMPH_EXCEPTION_LOCATION);
1452  }
1453  }
1454 
1455 
1456  //========================================================================
1457  /// Overload the derivative of the residuals and jacobian
1458  /// with respect to a parameter so that it breaks because it should not
1459  /// be required
1460  //========================================================================
1462  double* const &parameter_pt,
1463  Vector<double> &dres_dparam,
1464  DenseMatrix<double> &djac_dparam)
1465  {
1466  std::ostringstream error_stream;
1467  error_stream <<
1468  "This function has not been implemented because it is not required\n";
1469  error_stream << "in standard problems.\n";
1470  error_stream <<
1471  "If you find that you need it, you will have to implement it!\n\n";
1472 
1473  throw OomphLibError(error_stream.str(),
1474  OOMPH_CURRENT_FUNCTION,
1475  OOMPH_EXCEPTION_LOCATION);
1476  }
1477 
1478 
1479  //=====================================================================
1480  /// Overload the hessian vector product function so that
1481  /// it breaks because it should not be required
1482  //========================================================================
1484  GeneralisedElement* const &elem_pt,
1485  Vector<double> const &Y,
1486  DenseMatrix<double> const &C,
1487  DenseMatrix<double> &product)
1488  {
1489  std::ostringstream error_stream;
1490  error_stream <<
1491  "This function has not been implemented because it is not required\n";
1492  error_stream << "in standard problems.\n";
1493  error_stream <<
1494  "If you find that you need it, you will have to implement it!\n\n";
1495 
1496  throw OomphLibError(error_stream.str(),
1497  OOMPH_CURRENT_FUNCTION,
1498  OOMPH_EXCEPTION_LOCATION);
1499  }
1500 
1501 
1502  //==========================================================================
1503  /// Return the eigenfunction(s) associated with the bifurcation that
1504  /// has been detected in bifurcation tracking problems
1505  //==========================================================================
1507  {
1508  //There is only one (real) null vector
1509  eigenfunction.resize(1);
1511  //Rebuild the vector
1512  eigenfunction[0].build(&dist,0.0);
1513  //Set the value to be the null vector
1514  for(unsigned n=0;n<Ndof;n++)
1515  {
1516  eigenfunction[0][n] = Y[n];
1517  }
1518  }
1519 
1520 
1521 
1522  //=======================================================================
1523  /// Destructor return the problem to its original (non-augmented) state
1524  //=======================================================================
1526  {
1527  //If we are using the block solver reset the problem's linear solver
1528  //to the original one
1529  AugmentedBlockFoldLinearSolver *block_fold_solver_pt
1530  = dynamic_cast<AugmentedBlockFoldLinearSolver*>(
1532 
1533  if(block_fold_solver_pt)
1534  {
1535  //Reset the problem's linear solver
1536  Problem_pt->linear_solver_pt() = block_fold_solver_pt->linear_solver_pt();
1537  //Delete the block solver
1538  delete block_fold_solver_pt;
1539  }
1540 
1541  //Resize the number of dofs
1542  Problem_pt->Dof_pt.resize(Ndof);
1544  Ndof,false);
1545  //Remove all previous sparse storage used during Jacobian assembly
1547  }
1548 
1549  //====================================================================
1550  /// Set to solve the augmented-by-one block system.
1551  //===================================================================
1553  {
1554  //Only bother to do anything if we haven't already set the flag
1556  {
1557  //If we were solving the system with the original jacobian add the
1558  //parameter
1560  {
1561  Problem_pt->Dof_pt.push_back(Parameter_pt);
1562  }
1563 
1564  //Restrict the problem to the standard variables and
1565  //the bifurcation parameter only
1566  Problem_pt->Dof_pt.resize(Ndof+1);
1568  Ndof+1,false);
1569  //Remove all previous sparse storage used during Jacobian assembly
1571  //Set the solve flag
1573  }
1574  }
1575 
1576 
1577  //====================================================================
1578  /// Set to solve the block system. The boolean flag specifies
1579  /// whether the block decomposition should use exactly the same jacobian
1580  //===================================================================
1582  {
1583  //Only bother to do anything if we haven't already set the flag
1585  {
1586  //Restrict the problem to the standard variables
1587  Problem_pt->Dof_pt.resize(Ndof);
1589  Ndof,false);
1590  //Remove all previous sparse storage used during Jacobian assembly
1592 
1593  //Set the solve flag
1595  }
1596  }
1597 
1598  //=================================================================
1599  /// Set to Solve non-block system
1600  //=================================================================
1602  {
1603  //Only do something if we are not solving the full system
1605  {
1606  //If we were solving the problem without any augmentation,
1607  //add the parameter again
1609  {Problem_pt->Dof_pt.push_back(Parameter_pt);}
1610 
1611  //Always add the additional unknowns back into the problem
1612  for(unsigned n=0;n<Ndof;n++)
1613  {
1614  Problem_pt->Dof_pt.push_back(&Y[n]);
1615  }
1616 
1617  // update the Dof distribution pt
1619  Ndof*2+1,false);
1620  //Remove all previous sparse storage used during Jacobian assembly
1622 
1624  }
1625  }
1626 
1627 
1628 
1629  //======================================================================
1630  /// Clean up the memory that may have been allocated by the solver
1631  //=====================================================================
1633  {
1634  if(B_pt!=0) {delete B_pt;}
1635  if(C_pt!=0) {delete C_pt;}
1636  if(D_pt!=0) {delete D_pt;}
1637  if(dJy_dparam_pt!=0) {delete dJy_dparam_pt;}
1638  }
1639 
1640  //===================================================================
1641  /// Use a block factorisation to solve the augmented system
1642  /// associated with a PitchFork bifurcation.
1643  //===================================================================
1645  DoubleVector &result)
1646  {
1647  std::cout << "Block pitchfork solve" << std::endl;
1648  //Locally cache the pointer to the handler.
1649  PitchForkHandler* handler_pt =
1650  static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());
1651 
1652  //Get the augmented distribution from the handler and use it as
1653  //the distribution of the linear solver
1654  LinearAlgebraDistribution aug_dist =
1655  handler_pt->Augmented_dof_distribution_pt;
1656  this->build_distribution(aug_dist);
1657 
1658  // If the result vector is not setup then die
1659  if(!result.built())
1660  {
1661  throw OomphLibError("Result vector must be built\n",
1662  OOMPH_CURRENT_FUNCTION,
1663  OOMPH_EXCEPTION_LOCATION);
1664  }
1665 
1666  //Store the distribution of the result, which is probably not in
1667  //the natural distribution of the augmented solver
1668  LinearAlgebraDistribution result_dist(result.distribution_pt());
1669 
1670  //Locally cache a pointer to the parameter
1671  double* const parameter_pt = handler_pt->Parameter_pt;
1672 
1673  //Firstly get the derivatives of the full augmented system with
1674  //respect to the parameter
1675 
1676  //Allocate storage for the derivatives of the residuals with respect
1677  //to the global parameter using the augmented distribution
1678  DoubleVector dRdparam;
1679  //Then get the appropriate derivatives which will come back in
1680  //some distribution
1681  problem_pt->get_derivative_wrt_global_parameter(parameter_pt,dRdparam);
1682  //Redistribute into the augmented distribution
1683  dRdparam.redistribute(&aug_dist);
1684 
1685  //Switch the handler to "block solver" mode (sort out distribution)
1686  handler_pt->solve_block_system();
1687 
1688  //Temporary vector for the result (I shouldn't have to set this up)
1689  DoubleVector x1;
1690 
1691  //We are going to do resolves using the underlying linear solver
1693  //Solve the first (standard) system Jx1 = R
1694  Linear_solver_pt->solve(problem_pt,x1);
1695 
1696  //Allocate storage for B, C and D which can be used in the resolve
1697  //and the derivatives of the jacobian/eigenvector product with
1698  //respect to the parameter
1699  if(B_pt!=0) {delete B_pt;}
1701  if(C_pt!=0) {delete C_pt;}
1703  if(D_pt!=0) {delete D_pt;}
1705  //Need this to be in the distribution of the dofs
1706  if(dJy_dparam_pt!=0) {delete dJy_dparam_pt;}
1707  dJy_dparam_pt = new DoubleVector(handler_pt->Dof_distribution_pt,0.0);
1708 
1709  //Get the symmetry vector from the handler should have the
1710  //Dof distribution
1711  DoubleVector psi = handler_pt->Psi;
1712 
1713  //Temporary vector for the rhs that is dR/dparam (augmented distribution)
1714  //DoubleVector F(Linear_solver_pt->distribution_pt(),0.0);
1715  DoubleVector F(handler_pt->Dof_distribution_pt);
1716  const unsigned n_dof_local = F.nrow_local();
1717 
1718  //f.nrow local copied from dRdparam
1719  for(unsigned n=0;n<n_dof_local;n++) {F[n] = dRdparam[n];}
1720  //Fill in the rhs that is dJy/dparam //dJy_dparam nrow local
1721 
1722  //In standard cases the offset here will be one
1723  unsigned offset=1;
1724 #ifdef OOMPH_HAS_MPI
1725  //If we are distributed and not on the first processor
1726  //then there is no offset and the eigenfunction is
1727  //directly after the standard dofs
1728  int my_rank = problem_pt->communicator_pt()->my_rank();
1729  if((problem_pt->distributed()) && (my_rank!=0)) {offset=0;}
1730 #endif
1731  for(unsigned n=0;n<n_dof_local;n++)
1732  {(*dJy_dparam_pt)[n] = dRdparam[n_dof_local + offset + n];}
1733 
1734  //Now resolve to find c and d
1735  //First, redistribute F and psi into the linear algebra distribution
1738 
1740  Linear_solver_pt->resolve(psi,*D_pt);
1741 
1742  //We can now construct various dot products
1743  double psi_d = psi.dot(*D_pt);
1744  double psi_c = psi.dot(*C_pt);
1745  double psi_x1 = psi.dot(x1);
1746 
1747  //Calculate another intermediate constant
1748  const double Psi = psi_d/psi_c;
1749 
1750  //Construct the second intermediate value,
1751  //assumption is that result has been
1752  //set to the current value of the residuals
1753  //First, redistribute into the Natural distribution of the augmented system
1754  result.redistribute(&aug_dist);
1755  //The required parameter is that corresponding to the dof, which
1756  //is only stored on the root processor
1757  double parameter_residual = result[n_dof_local];
1758 #ifdef OOMPH_HAS_MPI
1759  //Broadcast to all others, if we have a distributed problem
1760  if(problem_pt->distributed())
1761  {
1762  MPI_Bcast(&parameter_residual,1,MPI_DOUBLE,0,
1763  problem_pt->communicator_pt()->mpi_comm());
1764  }
1765 #endif
1766  //Now construct the value
1767  double x2 = (psi_x1 - parameter_residual)/psi_c;
1768 
1769  //Now construct the vectors that multiply the jacobian terms
1771  D_and_X1[0].build(Linear_solver_pt->distribution_pt(),0.0);
1772  D_and_X1[1].build(Linear_solver_pt->distribution_pt(),0.0);
1773  //Fill in the appropriate terms
1774  //Get the number of local dofs from the Linear_solver_pt distribution
1775  const unsigned n_dof_local_linear_solver =
1777  for(unsigned n=0;n<n_dof_local_linear_solver;n++)
1778  {
1779  const double C_ = (*C_pt)[n];
1780  D_and_X1[0][n] = (*D_pt)[n] - Psi*C_;
1781  D_and_X1[1][n] = x1[n] - x2*C_;
1782  }
1783 
1784  //Local storage for the result of the product terms
1785  Vector<DoubleVectorWithHaloEntries> Jprod_D_and_X1(2);
1786 
1787  //Redistribute the inputs to have the Dof distribution pt
1788  D_and_X1[0].redistribute(handler_pt->Dof_distribution_pt);
1789  D_and_X1[1].redistribute(handler_pt->Dof_distribution_pt);
1790 
1791  //Set up the halo scheme
1792 #ifdef OOMPH_HAS_MPI
1793  D_and_X1[0].build_halo_scheme(problem_pt->Halo_scheme_pt);
1794  D_and_X1[1].build_halo_scheme(problem_pt->Halo_scheme_pt);
1795 #endif
1796 
1797  //Get the products from the new problem function
1798  problem_pt->get_hessian_vector_products(
1799  handler_pt->Y,D_and_X1,Jprod_D_and_X1);
1800 
1801  //OK, now we can formulate the next vectors
1802  //(again assuming result contains residuals)
1803  //Need to redistribute F to the dof distribution
1804  //F.redistribute(handler_pt->Dof_distribution_pt);
1805  DoubleVector G(handler_pt->Dof_distribution_pt);
1806 
1807  for(unsigned n=0;n<n_dof_local;n++)
1808  {
1809  G[n] = result[n_dof_local+offset+n]
1810  - Jprod_D_and_X1[1][n] - x2*dRdparam[n_dof_local+offset+n];
1811  Jprod_D_and_X1[0][n] *= -1.0;
1812  Jprod_D_and_X1[0][n] -= Psi*dRdparam[n_dof_local+offset+n];
1813  }
1814 
1815 
1816  //Then redistribute back to the linear solver distribution
1818  Jprod_D_and_X1[0].redistribute(Linear_solver_pt->distribution_pt());
1819  Jprod_D_and_X1[1].redistribute(Linear_solver_pt->distribution_pt());
1820 
1821  //Linear solve to get B
1822  Linear_solver_pt->resolve(Jprod_D_and_X1[0],*B_pt);
1823  //Liner solve to get x3
1825  Linear_solver_pt->resolve(G,x3);
1826 
1827  //Construst a couple of additional products
1828  double l_x3 = psi.dot(x3);
1829  double l_b = psi.dot(*B_pt);
1830 
1831  //get the last intermediate variable
1832  //The required parameter is that corresponding to the dof, which
1833  //is only stored on the root processor
1834  double sigma_residual = result[2*(n_dof_local+offset)-1];
1835 #ifdef OOMPH_HAS_MPI
1836  //Broadcast to all others, if we have a distributed problem
1837  if(problem_pt->distributed())
1838  {
1839  MPI_Bcast(&sigma_residual,1,MPI_DOUBLE,0,
1840  problem_pt->communicator_pt()->mpi_comm());
1841  }
1842 #endif
1843 
1844 
1845  const double delta_sigma = (l_x3 - sigma_residual)/l_b;
1846  const double delta_lambda = x2 - delta_sigma*Psi;
1847 
1848  //Need to do some more rearrangements here because result is global
1849  //but the other vectors are not!
1850 
1851  //Create temporary DoubleVectors to hold the results
1854 
1855  for(unsigned n=0;n<n_dof_local_linear_solver;n++)
1856  {
1857  res1[n] =
1858  x1[n] - delta_lambda*(*C_pt)[n] - delta_sigma*(*D_pt)[n];
1859  res2[n] = x3[n] - delta_sigma*(*B_pt)[n];
1860  }
1861 
1862  //Now redistribute these into the Dof distribution pointer
1863  res1.redistribute(handler_pt->Dof_distribution_pt);
1864  res2.redistribute(handler_pt->Dof_distribution_pt);
1865 
1866  //Now can copy over into results into the result vector
1867  for(unsigned n=0;n<n_dof_local;n++)
1868  {
1869  result[n] = res1[n];
1870  result[n_dof_local + offset + n] = res2[n];
1871  }
1872 
1873  //Add the final contributions to the residuals
1874  //only on the root processor if we have a distributed problem
1875 #ifdef OOMPH_HAS_MPI
1876  if((!problem_pt->distributed()) || (my_rank==0))
1877 #endif
1878  {
1879  result[n_dof_local] = delta_lambda;
1880  result[2*n_dof_local+1] = delta_sigma;
1881  }
1882 
1883 
1884  //The sign of the determinant is given by the sign of
1885  //the product psi_c and l_b
1886  //NOT CHECKED YET!
1887  problem_pt->sign_of_jacobian() =
1888  static_cast<int>(std::fabs(psi_c*l_b)/(psi_c*l_b));
1889 
1890  //Redistribute the result into its incoming distribution
1891  result.redistribute(&result_dist);
1892 
1893  //Switch things to our block solver
1894  handler_pt->solve_full_system();
1895 
1896  //If we are not storing things, clear the results
1897  if(!Enable_resolve)
1898  {
1899  //We no longer need to store the matrix
1901  delete B_pt; B_pt=0;
1902  delete C_pt; C_pt=0;
1903  delete D_pt; D_pt=0;
1904  delete dJy_dparam_pt; dJy_dparam_pt = 0;
1905  }
1906  //Otherwise also store the pointer to the problem
1907  else
1908  {
1909  Problem_pt = problem_pt;
1910  }
1911  }
1912 
1913 
1914  //==============================================================
1915  //Hack the re-solve to use the block factorisation
1916  //==============================================================
1918  DoubleVector &result)
1919  {
1920  std::cout << "Block pitchfork resolve" << std::endl;
1921  //Check that the factors have been stored
1922  if(B_pt==0)
1923  {
1924  throw OomphLibError("The required vectors have not been stored",
1925  OOMPH_CURRENT_FUNCTION,
1926  OOMPH_EXCEPTION_LOCATION);
1927  }
1928 
1929 
1930  //Cache pointer to the problem
1931  Problem* const problem_pt = Problem_pt;
1932 
1933  //Locally cache pointer to the handler
1934  PitchForkHandler* handler_pt =
1935  static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());
1936 
1937  //Get the augmented distribution from the handler
1938  LinearAlgebraDistribution aug_dist =
1939  handler_pt->Augmented_dof_distribution_pt;
1940  //this->build_distribution(aug_dist);
1941 
1942  //Find the number of dofs of the augmented system
1943  //const unsigned n_aug_dof = problem_pt->ndof();
1944 
1945  //Storage for the result distribution
1946  LinearAlgebraDistribution result_dist;
1947 
1948  // if the result vector is not setup then rebuild with distribution
1949  // = natural distribution of augmented solver
1950  if(!result.built())
1951  {
1952  result.build(&aug_dist,0.0);
1953  }
1954  //Otherwise store the incoming distribution and redistribute
1955  else
1956  {
1957  result_dist.build(result.distribution_pt());
1958  result.redistribute(&aug_dist);
1959  }
1960 
1961 
1962  //Locally cache a pointer to the parameter
1963  //double* const parameter_pt = handler_pt->Parameter_pt;
1964 
1965  //Switch things to our block solver
1966  handler_pt->solve_block_system();
1967  //We need to find out the number of dofs
1968  //unsigned n_dof = problem_pt->ndof();
1969 
1970  // create the linear algebra distribution for this solver
1971  // currently only global (non-distributed) distributions are allowed
1972  //LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
1973  //this->build_distribution(dist);
1974 
1975  // if the result vector is not setup then rebuild with distribution = global
1976  //if (!result.built())
1977  // {
1978  // result.build(this->distribution_pt(),0.0);
1979  // }
1980 
1981 
1982  //Copy the rhs into local storage
1983  DoubleVector rhs_local = rhs;
1984  //and redistribute into the augmented distribution
1985  rhs_local.redistribute(&aug_dist);
1986 
1987  //Setup storage for output
1990 
1991  //Create input for RHS with the natural distribution
1992  DoubleVector input_x1(handler_pt->Dof_distribution_pt);
1993  const unsigned n_dof_local = input_x1.nrow_local();
1994 
1995  //Set the values of the a vector
1996  for(unsigned n=0;n<n_dof_local;n++) {input_x1[n] = rhs_local[n];}
1997  //Need to redistribute into the linear algebra distribution
1999 
2000  //We want to resolve, of course
2002  //Now solve for the first vector
2003  Linear_solver_pt->resolve(input_x1,x1);
2004 
2005  //Get the symmetry vector from the handler
2006  DoubleVector psi = handler_pt->Psi;
2007  //redistribute local copy
2009 
2010  //We can now construct various dot products
2011  double psi_d = psi.dot(*D_pt);
2012  double psi_c = psi.dot(*C_pt);
2013  double psi_x1 = psi.dot(x1);
2014 
2015  //Calculate another intermediate constant
2016  const double Psi = psi_d/psi_c;
2017 
2018  //Construct the second intermediate value,
2019  //assumption is that rhs has been set to the current value of the residuals
2020  double parameter_residual = rhs_local[n_dof_local];
2021 #ifdef OOMPH_HAS_MPI
2022  //Broadcast to all others, if we have a distributed problem
2023  if(problem_pt->distributed())
2024  {
2025  MPI_Bcast(&parameter_residual,1,MPI_DOUBLE,0,
2026  problem_pt->communicator_pt()->mpi_comm());
2027  }
2028 #endif
2029 
2030  double x2 = (psi_x1 - parameter_residual)/psi_c;
2031 
2032  //Now construct the vectors that multiply the jacobian terms
2033  //Vector<double> X1(n_dof/*+1*/);
2035  X1[0].build(Linear_solver_pt->distribution_pt(),0.0);
2036 
2037  const unsigned n_dof_local_linear_solver =
2039  for(unsigned n=0;n<n_dof_local_linear_solver;n++)
2040  {
2041  X1[0][n] = x1[n] - x2*(*C_pt)[n];
2042  }
2043 
2044  //Local storage for the product term
2046 
2047  //Redistribute the inputs to have the Dof distribution pt
2048  X1[0].redistribute(handler_pt->Dof_distribution_pt);
2049 
2050  //Set up the halo scheme
2051 #ifdef OOMPH_HAS_MPI
2052  X1[0].build_halo_scheme(problem_pt->Halo_scheme_pt);
2053 #endif
2054 
2055  //Get the product from the problem
2056  problem_pt->get_hessian_vector_products(handler_pt->Y,X1,Jprod_X1);
2057 
2058  //In standard cases the offset here will be one
2059  unsigned offset=1;
2060 #ifdef OOMPH_HAS_MPI
2061  //If we are distributed and not on the first processor
2062  //then there is no offset and the eigenfunction is
2063  //directly after the standard dofs
2064  int my_rank = problem_pt->communicator_pt()->my_rank();
2065  if((problem_pt->distributed()) && (my_rank!=0)) {offset=0;}
2066 #endif
2067 
2068  //OK, now we can formulate the next vectors
2069  //(again assuming result contains residuals)
2070  //Local storage for the product terms
2071  DoubleVector Mod_Jprod_X1(handler_pt->Dof_distribution_pt,0.0);
2072 
2073  for(unsigned n=0;n<n_dof_local;n++)
2074  {
2075  Mod_Jprod_X1[n] = rhs_local[n_dof_local+offset+n] - Jprod_X1[0][n]
2076  - x2*(*dJy_dparam_pt)[n];
2077  }
2078 
2079  //Redistribute back to the linear solver distribution
2081 
2082  //Liner solve to get x3
2083  Linear_solver_pt->resolve(Mod_Jprod_X1,x3);
2084 
2085  //Construst a couple of additional products
2086  double l_x3 = psi.dot(x3);
2087  double l_b = psi.dot(*B_pt);
2088 
2089  //get the last intermediate variable
2090  //The required parameter is that corresponding to the dof, which
2091  //is only stored on the root processor
2092  double sigma_residual = rhs_local[2*(n_dof_local+offset)-1];
2093 #ifdef OOMPH_HAS_MPI
2094  //Broadcast to all others, if we have a distributed problem
2095  if(problem_pt->distributed())
2096  {
2097  MPI_Bcast(&sigma_residual,1,MPI_DOUBLE,0,
2098  problem_pt->communicator_pt()->mpi_comm());
2099  }
2100 #endif
2101 
2102  //get the last intermediate variable
2103  const double delta_sigma = (l_x3 - sigma_residual)/l_b;
2104  const double delta_lambda = x2 - delta_sigma*Psi;
2105 
2106  //Create temporary DoubleVectors to hold the results
2109 
2110  for(unsigned n=0;n<n_dof_local_linear_solver;n++)
2111  {
2112  res1[n] = x1[n] - delta_lambda*(*C_pt)[n] - delta_sigma*(*D_pt)[n];
2113  res2[n] = x3[n] - delta_sigma*(*B_pt)[n];
2114  }
2115 
2116  //Now redistribute these into the Dof distribution pointer
2117  res1.redistribute(handler_pt->Dof_distribution_pt);
2118  res2.redistribute(handler_pt->Dof_distribution_pt);
2119 
2120  //Now can copy over into results into the result vector
2121  for(unsigned n=0;n<n_dof_local;n++)
2122  {
2123  result[n] = res1[n];
2124  result[n_dof_local + offset + n] = res2[n];
2125  }
2126 
2127  //Add the final contributions to the residuals
2128  //only on the root processor if we have a distributed problem
2129 #ifdef OOMPH_HAS_MPI
2130  if((!problem_pt->distributed()) || (my_rank==0))
2131 #endif
2132  {
2133  result[n_dof_local] = delta_lambda;
2134  result[2*n_dof_local+1] = delta_sigma;
2135  }
2136 
2137  //Redistribute the result into its incoming distribution (if it had one)
2138  if(result_dist.built()) {result.redistribute(&result_dist);}
2139 
2141 
2142  //Switch things to our block solver
2143  handler_pt->solve_full_system();
2144  }
2145 
2146 
2147 //--------------------------------------------------------------
2148 
2149 
2150 
2151  //======================================================================
2152  /// Clean up the memory that may have been allocated by the solver
2153  //=====================================================================
2155  {
2156  if(Alpha_pt!=0) {delete Alpha_pt;}
2157  if(E_pt!=0) {delete E_pt;}
2158  }
2159 
2160  //===================================================================
2161  /// Use a block factorisation to solve the augmented system
2162  /// associated with a PitchFork bifurcation.
2163  //===================================================================
2165  DoubleVector &result)
2166  {
2167  std::cout << "Augmented pitchfork solve" << std::endl;
2168  // if the result is setup then it should not be distributed
2169 #ifdef PARANOID
2170  if (result.built())
2171  {
2172  if (result.distributed())
2173  {
2174  throw OomphLibError("The result vector must not be distributed",
2175  OOMPH_CURRENT_FUNCTION,
2176  OOMPH_EXCEPTION_LOCATION);
2177  }
2178  }
2179 #endif
2180 
2181 
2182  //Locally cache the pointer to the handler.
2183  PitchForkHandler* handler_pt =
2184  static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());
2185 
2186  //Switch the handler to "block solver" mode
2187  handler_pt->solve_augmented_block_system();
2188 
2189  //We need to find out the number of dofs in the problem
2190  unsigned n_dof = problem_pt->ndof();
2191 
2192  // create the linear algebra distribution for this solver
2193  // currently only global (non-distributed) distributions are allowed
2194  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
2195  this->build_distribution(dist);
2196 
2197  // if the result vector is not setup then rebuild with distribution = global
2198  if (!result.built())
2199  {
2200  result.build(this->distribution_pt(),0.0);
2201  }
2202 
2203  //Setup storage for temporary vectors
2204  DoubleVector a(this->distribution_pt(),0.0), b(this->distribution_pt(),0.0);
2205 
2206  //Allocate storage for Alpha which can be used in the resolve
2207  if(Alpha_pt!=0) {delete Alpha_pt;}
2208  Alpha_pt = new DoubleVector(this->distribution_pt(),0.0);
2209 
2210  //We are going to do resolves using the underlying linear solver
2212  //Solve the first system Aa = R
2213  Linear_solver_pt->solve(problem_pt,a);
2214 
2215  //Get the symmetry vector from the handler
2216  DoubleVector psi(this->distribution_pt(),0.0);
2217  for(unsigned n=0;n<(n_dof-1);++n) {psi[n] = handler_pt->Psi[n];}
2218  //Set the final entry to zero
2219  psi[n_dof-1] = 0.0;
2220 
2221  //Now resolve to find alpha
2223 
2224  //We can now construct our multipliers
2225  //Prepare to scale
2226  double dof_length=0.0, a_length=0.0, alpha_length=0.0;
2227  for(unsigned n=0;n<n_dof;n++)
2228  {
2229  if(std::fabs(problem_pt->dof(n)) > dof_length)
2230  {dof_length = std::fabs(problem_pt->dof(n));}
2231  if(std::fabs(a[n]) > a_length) {a_length = std::fabs(a[n]);}
2232  if(std::fabs((*Alpha_pt)[n]) > alpha_length)
2233  {alpha_length = std::fabs((*Alpha_pt)[n]);}
2234  }
2235 
2236  double a_mult = dof_length/a_length;
2237  double alpha_mult = dof_length/alpha_length;
2238  const double FD_step = 1.0e-8;
2239  a_mult += FD_step; alpha_mult += FD_step;
2240  a_mult *= FD_step; alpha_mult *= FD_step;
2241 
2242  //Local storage for the product terms
2243  DoubleVector Jprod_a(this->distribution_pt(),0.0),
2244  Jprod_alpha(this->distribution_pt(),0.0);
2245 
2246  //Calculate the product of the jacobian matrices, etc
2247  unsigned long n_element = problem_pt->mesh_pt()->nelement();
2248  for(unsigned long e = 0;e<n_element;e++)
2249  {
2250  GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
2251  //Loop over the ndofs in each element
2252  unsigned n_var = handler_pt->ndof(elem_pt);
2253  //Get the jacobian matrices
2254  DenseMatrix<double> jac(n_var), jac_a(n_var), jac_alpha(n_var);
2255  // the elemental residual
2256  Vector<double> res_elemental(n_var);
2257  //Get unperturbed jacobian
2258  handler_pt->get_jacobian(elem_pt,res_elemental,jac);
2259 
2260  //Backup the dofs
2261  Vector<double> dof_bac(n_var);
2262  //Perturb the dofs
2263  for(unsigned n=0;n<n_var;n++)
2264  {
2265  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
2266  dof_bac[n] = problem_pt->dof(eqn_number);
2267  //Pertub by vector a
2268  problem_pt->dof(eqn_number) += a_mult*a[eqn_number];
2269  }
2270 
2272 
2273  //Now get the new jacobian
2274  handler_pt->get_jacobian(elem_pt,res_elemental,jac_a);
2275 
2276  //Perturb the dofs
2277  for(unsigned n=0;n<n_var;n++)
2278  {
2279  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
2280  problem_pt->dof(eqn_number) = dof_bac[n];
2281  //Pertub by vector a
2282  problem_pt->dof(eqn_number) += alpha_mult*(*Alpha_pt)[eqn_number];
2283  }
2284 
2286 
2287  //Now get the new jacobian
2288  handler_pt->get_jacobian(elem_pt,res_elemental,jac_alpha);
2289 
2290  //Reset the dofs
2291  for(unsigned n=0;n<n_var;n++)
2292  {
2293  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
2294  problem_pt->dof(eqn_number) = dof_bac[n];
2295  }
2296 
2298 
2299  //OK, now work out the products
2300  for(unsigned n=0;n<(n_var-1);n++)
2301  {
2302  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
2303  double prod_a=0.0, prod_alpha=0.0;
2304  for(unsigned m=0;m<(n_var-1);m++)
2305  {
2306  unsigned unknown = handler_pt->eqn_number(elem_pt,m);
2307  prod_a += (jac_a(n,m) - jac(n,m))*
2308  handler_pt->Y[unknown];
2309  prod_alpha += (jac_alpha(n,m) - jac(n,m))*
2310  handler_pt->Y[unknown];
2311  }
2312  Jprod_a[eqn_number] += prod_a/a_mult;
2313  Jprod_alpha[eqn_number] += prod_alpha/alpha_mult;
2314  }
2315  }
2316 
2317  Jprod_alpha[n_dof-1] = 0.0;
2318  Jprod_a[n_dof-1] = 0.0;
2319 
2320  //OK, now we can formulate the next vectors
2321  //The assumption here is that the result has been set to the
2322  //residuals.
2323  for(unsigned n=0;n<n_dof-1;n++)
2324  {
2325  b[n] = result[n_dof+n] - Jprod_a[n];
2326  }
2327  b[n_dof-1] = result[2*n_dof-1];
2328 
2329 
2330  //Allocate storage for E which can be used in the resolve
2331  if(E_pt!=0) {delete E_pt;}
2332  E_pt = new DoubleVector(this->distribution_pt(),0.0);
2333 
2334  DoubleVector f(this->distribution_pt(),0.0);
2335 
2336  Linear_solver_pt->resolve(b,f);
2337  Linear_solver_pt->resolve(Jprod_alpha,*E_pt);
2338 
2339  //Calculate the final entry in the vector e
2340  const double e_final = (*E_pt)[n_dof-1];
2341  //Calculate the final entry in the vector d
2342  const double d_final = -f[n_dof-1]/e_final;
2343  //Assemble the final corrections
2344  for(unsigned n=0;n<n_dof-1;n++)
2345  {
2346  result[n] = a[n] - (*Alpha_pt)[n]*d_final;
2347  result[n_dof+n] = f[n] + (*E_pt)[n]*d_final;
2348  }
2349 
2350  result[n_dof-1] = a[n_dof-1] - (*Alpha_pt)[n_dof-1]*d_final;
2351  result[2*n_dof-1] = d_final;
2352 
2353  //The sign of the jacobian is the sign of the final entry in e
2354  problem_pt->sign_of_jacobian() =
2355  static_cast<int>(std::fabs(e_final)/e_final);
2356 
2357 
2358  //Switch things to our block solver
2359  handler_pt->solve_full_system();
2360 
2361  //If we are not storing things, clear the results
2362  if(!Enable_resolve)
2363  {
2364  //We no longer need to store the matrix
2366  delete Alpha_pt; Alpha_pt=0;
2367  delete E_pt; E_pt=0;
2368  }
2369  //Otherwise also store the pointer to the problem
2370  else
2371  {
2372  Problem_pt = problem_pt;
2373  }
2374  }
2375 
2376 
2377  //==============================================================
2378  //Hack the re-solve to use the block factorisation
2379  //==============================================================
2381  DoubleVector &result)
2382  {
2383  std::cout << "Augmented pitchfork resolve" << std::endl;
2384  //Check that the factors have been stored
2385  if(Alpha_pt==0)
2386  {
2387  throw OomphLibError("The required vectors have not been stored",
2388  OOMPH_CURRENT_FUNCTION,
2389  OOMPH_EXCEPTION_LOCATION);
2390  }
2391 
2392  //Get the pointer to the problem
2393  Problem* const problem_pt = Problem_pt;
2394 
2395  PitchForkHandler* handler_pt =
2396  static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());
2397 
2398  //Switch things to our block solver
2399  handler_pt->solve_augmented_block_system();
2400  //We need to find out the number of dofs
2401  unsigned n_dof = problem_pt->ndof();
2402 
2403  // create the linear algebra distribution for this solver
2404  // currently only global (non-distributed) distributions are allowed
2405  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
2406  this->build_distribution(dist);
2407 
2408  // if the result vector is not setup then rebuild with distribution = global
2409  if (!result.built())
2410  {
2411  result.build(this->distribution_pt(),0.0);
2412  }
2413 
2414 
2415  //Setup storage
2416  DoubleVector a(this->distribution_pt(),0.0), b(this->distribution_pt(),0.0);
2417 
2418  //Set the values of the a vector
2419  for(unsigned n=0;n<n_dof;n++) {a[n] = rhs[n];}
2420 
2422 
2423  // Copy rhs vector into local storage so it doesn't get overwritten
2424  // if the linear solver decides to initialise the solution vector, say,
2425  // which it's quite entitled to do!
2426  DoubleVector input_a(a);
2427 
2428  Linear_solver_pt->resolve(input_a,a);
2429 
2430  //We can now construct our multipliers
2431  //Prepare to scale
2432  double dof_length=0.0, a_length=0.0;
2433  for(unsigned n=0;n<n_dof;n++)
2434  {
2435  if(std::fabs(problem_pt->dof(n)) > dof_length)
2436  {dof_length = std::fabs(problem_pt->dof(n));}
2437 
2438  if(std::fabs(a[n]) > a_length) {a_length = std::fabs(a[n]);}
2439  }
2440  double a_mult = dof_length/a_length;
2441  const double FD_step = 1.0e-8;
2442  a_mult += FD_step;
2443  a_mult *= FD_step;
2444 
2445  DoubleVector Jprod_a(this->distribution_pt(),0.0);
2446 
2447  unsigned long n_element = problem_pt->mesh_pt()->nelement();
2448  for(unsigned long e = 0;e<n_element;e++)
2449  {
2450  GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
2451  //Loop over the ndofs in each element
2452  unsigned n_var = handler_pt->ndof(elem_pt);
2453  //Get some jacobian matrices
2454  DenseMatrix<double> jac(n_var), jac_a(n_var);
2455  // the elemental residual
2456  Vector<double> res_elemental(n_var);
2457  //Get unperturbed jacobian
2458  handler_pt->get_jacobian(elem_pt,res_elemental,jac);
2459 
2460  //Backup the dofs
2461  DoubleVector dof_bac(this->distribution_pt(),0.0);
2462  //Perturb the dofs
2463  for(unsigned n=0;n<n_var;n++)
2464  {
2465  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
2466  dof_bac[n] = problem_pt->dof(eqn_number);
2467  //Pertub by vector a
2468  problem_pt->dof(eqn_number) += a_mult*a[eqn_number];
2469  }
2470 
2472 
2473  //Now get the new jacobian
2474  handler_pt->get_jacobian(elem_pt,res_elemental,jac_a);
2475 
2476  //Reset the dofs
2477  for(unsigned n=0;n<n_var;n++)
2478  {
2479  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
2480  problem_pt->dof(eqn_number) = dof_bac[n];
2481  }
2482 
2484 
2485  //OK, now work out the products
2486  for(unsigned n=0;n<(n_var-1);n++)
2487  {
2488  unsigned eqn_number = handler_pt->eqn_number(elem_pt,n);
2489  double prod_a=0.0;
2490  for(unsigned m=0;m<(n_var-1);m++)
2491  {
2492  unsigned unknown = handler_pt->eqn_number(elem_pt,m);
2493  prod_a += (jac_a(n,m) - jac(n,m))*
2494  handler_pt->Y[unknown];
2495  }
2496  Jprod_a[eqn_number] += prod_a/a_mult;
2497  }
2498  }
2499 
2500  Jprod_a[n_dof-1] = 0.0;
2501 
2502  //OK, now we can formulate the next vectors
2503  for(unsigned n=0;n<n_dof-1;n++)
2504  {
2505  b[n] = rhs[n_dof+n] - Jprod_a[n];
2506  }
2507  b[n_dof-1] = rhs[2*n_dof-1];
2508 
2509  DoubleVector f(this->distribution_pt(),0.0);
2510 
2511  Linear_solver_pt->resolve(b,f);
2512 
2513  //Calculate the final entry in the vector d
2514  const double d_final = -f[n_dof-1]/(*E_pt)[n_dof-1];
2515  //Assemble the final corrections
2516  for(unsigned n=0;n<n_dof-1;n++)
2517  {
2518  result[n] = a[n] - (*Alpha_pt)[n]*d_final;
2519  result[n_dof+n] = f[n] + (*E_pt)[n]*d_final;
2520  }
2521 
2522  result[n_dof-1] = a[n_dof-1] - (*Alpha_pt)[n_dof-1]*d_final;
2523  result[2*n_dof-1] = d_final;
2524 
2526 
2527  //Switch things to our block solver
2528  handler_pt->solve_full_system();
2529  }
2530 
2531  ///////////////////////////////////////////////////////////////////////
2532  // Non-inline functions for the PitchForkHandler class
2533  //////////////////////////////////////////////////////////////////////
2534 
2535  //==================================================================
2536  ///Constructor: Initialise the PitchForkHandler by setting intial
2537  ///guesses for Sigma, Y, specifying C and Psi and calculating count.
2538  //==================================================================
2540  AssemblyHandler* const
2541  &assembly_handler_pt,
2542  double* const &parameter_pt,
2543  const DoubleVector &symmetry_vector) :
2544  Solve_which_system(Full_augmented),
2545  Sigma(0.0),
2546  Parameter_pt(parameter_pt)
2547  {
2548  //Set the problem pointer
2549  Problem_pt = problem_pt;
2550  //Set the assembly handler
2551  Assembly_handler_pt = assembly_handler_pt;
2552  //Set the number of degrees of freedom
2553  Ndof = Problem_pt->ndof();
2554  //Backup the distribution
2556 #ifdef OOMPH_HAS_MPI
2557  //Set the distribution flag
2559 #endif
2560 
2561  //Find the elements in the problem
2562  unsigned n_element = Problem_pt->mesh_pt()->nelement();
2563 
2564 #ifdef OOMPH_HAS_MPI
2565 
2566  //Work out all the global equations to which this processor
2567  //contributes and store the halo data in the problem
2569 
2570 #endif
2571 
2572 
2573  //Now use the dof distribution for all double vectors
2578 #ifdef OOMPH_HAS_MPI
2580 #endif
2581 
2582  //Loop over the elements and count the entries
2583  //and number of (non-halo) elements
2584  unsigned n_non_halo_element_local=0;
2585  for(unsigned e=0;e<n_element;e++)
2586  {
2588 #ifdef OOMPH_HAS_MPI
2589  //Ignore halo elements
2590  if(!elem_pt->is_halo())
2591  {
2592 #endif
2593  //Increment the number of non halo elements
2594  ++n_non_halo_element_local;
2595  //Now count the number of times the element contributes to a value
2596  unsigned n_var = assembly_handler_pt->ndof(elem_pt);
2597  for(unsigned n=0;n<n_var;n++)
2598  {
2599  ++Count.global_value(assembly_handler_pt->eqn_number(elem_pt,n));
2600  }
2601 #ifdef OOMPH_HAS_MPI
2602  }
2603 #endif
2604  }
2605 
2606  //Add together all the counts
2607 #ifdef OOMPH_HAS_MPI
2609 
2610  //If distributed, find the total number of elements in the problem
2611  if(Distributed)
2612  {
2613  //Need to gather the total number of non halo elements
2614  MPI_Allreduce(&n_non_halo_element_local,&Nelement,1,MPI_UNSIGNED,MPI_SUM,
2615  Problem_pt->communicator_pt()->mpi_comm());
2616  }
2617  //Otherwise the total number is the same on each processor
2618  else
2619 #endif
2620  {
2621  Nelement = n_non_halo_element_local;
2622  }
2623 
2624 #ifdef OOMPH_HAS_MPI
2625  //Only add the parameter to the first processor, if distributed
2626  int my_rank = Problem_pt->communicator_pt()->my_rank();
2627  if((!Distributed) || (my_rank==0))
2628 #endif
2629  {
2630  //Add the parameter to the problem
2631  Problem_pt->Dof_pt.push_back(parameter_pt);
2632  }
2633 
2634  //Find length of the symmetry vector
2635  double length = symmetry_vector.norm();
2636 
2637  //Add the unknowns for the null vector to the problem
2638  //Normalise the symmetry vector and initialise the null vector to the
2639  //symmetry vector and set the constant vector, c, to the
2640  //normalised symmetry vector.
2641 
2642  //Need to redistribute the symmetry vector into the (natural)
2643  //distribution of the Dofs
2644 #ifdef OOMPH_HAS_MPI
2645  //Check that the symmetry vector has the correct distribution
2646  if(!(*symmetry_vector.distribution_pt() == *Dof_distribution_pt))
2647  {
2648  throw OomphLibError
2649  ("The symmetry vector must have the same distribution as the dofs\n",
2650  OOMPH_CURRENT_FUNCTION,
2651  OOMPH_EXCEPTION_LOCATION);
2652  }
2653 #endif
2654 
2655  //Only loop over the local unknowns
2656  const unsigned n_dof_local = Dof_distribution_pt->nrow_local();
2657  for(unsigned n=0;n<n_dof_local;n++)
2658  {
2659  Problem_pt->Dof_pt.push_back(&Y[n]);
2660  //Psi[n] = symmetry_vector[n];
2661  Psi[n] = Y[n] = C[n] = symmetry_vector[n]/length;
2662  }
2663 
2664 
2665 #ifdef OOMPH_HAS_MPI
2666  //Set up the required halo schemes (which also synchronises)
2670 
2671 
2672  if((!Distributed) || (my_rank==0))
2673 #endif
2674  //Add the slack parameter to the problem on the first processor
2675  //if distributed
2676  {
2677  Problem_pt->Dof_pt.push_back(&Sigma);
2678  }
2679 
2680 #ifdef OOMPH_HAS_MPI
2681  if(Distributed)
2682  {
2683  unsigned augmented_first_row=0;
2684  unsigned augmented_n_row_local=0;
2685 
2686  //Set up the translation scheme on every processor
2687  Global_eqn_number.resize(2*Ndof+2);
2688  int n_proc = Problem_pt->communicator_pt()->nproc();
2689  unsigned global_eqn_count=0;
2690  for(int d=0;d<n_proc;d++)
2691  {
2692  //Find out the first row of the current processor
2693  if(my_rank==d) {augmented_first_row = global_eqn_count;}
2694 
2695  const unsigned n_row_local = Dof_distribution_pt->nrow_local(d);
2696  const unsigned first_row = Dof_distribution_pt->first_row(d);
2697  //Add the basic equations
2698  for(unsigned n=0;n<n_row_local;n++)
2699  {
2700  Global_eqn_number[first_row+n] = global_eqn_count;
2701  ++global_eqn_count;
2702  }
2703  //If on the first processor add the pointer to the parameter
2704  if(d==0)
2705  {
2706  Global_eqn_number[Ndof] = global_eqn_count;
2707  ++global_eqn_count;
2708  }
2709  //Add the eigenfunction
2710  for(unsigned n=0;n<n_row_local;n++)
2711  {
2712  Global_eqn_number[Ndof + 1 + first_row+n] = global_eqn_count;
2713  ++global_eqn_count;
2714  }
2715  //Finally add the slack parameter
2716  if(d==0)
2717  {
2718  Global_eqn_number[2*Ndof+1] = global_eqn_count;
2719  ++global_eqn_count;
2720  }
2721  //Find out the number of rows of the current processor
2722  if(my_rank==d)
2723  {augmented_n_row_local = global_eqn_count - augmented_first_row;}
2724  }
2725 
2726  //Make a new linear algebra distribution
2729  augmented_first_row,
2730  augmented_n_row_local);
2731  }
2732  else
2733 #endif
2734  {
2737  2*Ndof+2,false);
2738  }
2739 
2740  // resize
2741  //Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
2742  // 2*Ndof+2,false);
2743 
2745 
2746  //Remove all previous sparse storage used during Jacobian assembly
2748  }
2749 
2750  //=====================================================================
2751  /// Destructor return the problem to its original (non-augmented) state
2752  //=====================================================================
2754  {
2755  //If we are using the block solver reset the problem's linear solver
2756  //to the original one
2757  BlockPitchForkLinearSolver* block_pitchfork_solver_pt
2758  = dynamic_cast<BlockPitchForkLinearSolver*>(
2760 
2761  if(block_pitchfork_solver_pt)
2762  {
2763  //Reset the problem's linear solver
2765  block_pitchfork_solver_pt->linear_solver_pt();
2766  //Delete the block solver
2767  delete block_pitchfork_solver_pt;
2768  }
2769 
2770  //If we are using the augmented
2771  //block solver reset the problem's linear solver
2772  //to the original one
2773  AugmentedBlockPitchForkLinearSolver* augmented_block_pitchfork_solver_pt
2774  = dynamic_cast<AugmentedBlockPitchForkLinearSolver*>(
2776 
2777  if(augmented_block_pitchfork_solver_pt)
2778  {
2779  //Reset the problem's linear solver
2781  augmented_block_pitchfork_solver_pt->linear_solver_pt();
2782  //Delete the block solver
2783  delete augmented_block_pitchfork_solver_pt;
2784  }
2785 
2786  //Now return the problem to its original size
2787  Problem_pt->Dof_pt.resize(this->Dof_distribution_pt->nrow_local());
2788 
2789  //Problem_pt->Dof_pt.resize(Ndof);
2790  //Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
2791  // Ndof,false);
2792 
2794  //Remove all previous sparse storage used during Jacobian assembly
2796  }
2797 
2798  //================================================================
2799  ///Get the number of elemental degrees of freedom
2800  //================================================================
2801  unsigned PitchForkHandler::ndof(GeneralisedElement* const &elem_pt)
2802  {
2803  unsigned raw_ndof = elem_pt->ndof();
2804 
2805  //Return different values depending on the type of block decomposition
2806  switch(Solve_which_system)
2807  {
2808  case Full_augmented:
2809  return (2*raw_ndof + 2);
2810  break;
2811 
2812  case Block_augmented_J:
2813  return (raw_ndof+1);
2814  break;
2815 
2816  case Block_J:
2817  return raw_ndof;
2818  break;
2819 
2820  default:
2821  std::ostringstream error_stream;
2822  error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
2823  << " not " << Solve_which_system << "\n";
2824  throw OomphLibError(error_stream.str(),
2825  OOMPH_CURRENT_FUNCTION,
2826  OOMPH_EXCEPTION_LOCATION);
2827  }
2828  }
2829 
2830  ///Get the global equation number of the local unknown
2831  unsigned long PitchForkHandler::eqn_number(GeneralisedElement* const &elem_pt,
2832  const unsigned &ieqn_local)
2833  {
2834  //Get the raw value
2835  unsigned raw_ndof = elem_pt->ndof();
2836  unsigned long global_eqn=0;
2837 
2838  //Which system are we solving
2839  switch(Solve_which_system)
2840  {
2841  //In the block case, it's just the standard numbering
2842  case Block_J:
2843  global_eqn = elem_pt->eqn_number(ieqn_local);
2844  break;
2845 
2846  //Block augmented not done properly yet
2847  case Block_augmented_J:
2848  //Not done the distributed case
2849 #ifdef OOMPH_HAS_MPI
2850  if(Distributed)
2851  {
2852  throw OomphLibError(
2853  "Block Augmented solver not implemented for distributed case\n",
2854  OOMPH_CURRENT_FUNCTION,
2855  OOMPH_EXCEPTION_LOCATION);
2856  }
2857 #endif
2858 
2859  //Full case
2860  case Full_augmented:
2861  //The usual equations
2862  if(ieqn_local < raw_ndof)
2863  {
2864  global_eqn = this->global_eqn_number(elem_pt->eqn_number(ieqn_local));
2865  }
2866  //The bifurcation parameter equation
2867  else if(ieqn_local == raw_ndof)
2868  {
2869  global_eqn = this->global_eqn_number(Ndof);
2870  }
2871  //If we are assembling the full system we also have
2872  //The components of the null vector
2873  else if(ieqn_local < (2*raw_ndof + 1))
2874  {
2875  global_eqn = this->global_eqn_number(
2876  Ndof + 1 + elem_pt->eqn_number(ieqn_local - 1 - raw_ndof));
2877  }
2878  //The slack parameter
2879  else
2880  {
2881  global_eqn = this->global_eqn_number(2*Ndof+1);
2882  }
2883  break;
2884  } //End of switch
2885  return global_eqn;
2886  }
2887 
2888  //==============================================================
2889  ///Get the residuals
2890  //==============================================================
2892  Vector<double> &residuals)
2893  {
2894 
2895  //Need to get raw residuals and jacobian
2896  unsigned raw_ndof = elem_pt->ndof();
2897 
2898  //Find out which system we are solving
2899  switch(Solve_which_system)
2900  {
2901  //If we are solving the original system
2902  case Block_J:
2903  {
2904  //get the basic residuals
2905  elem_pt->get_residuals(residuals);
2906  //Now multiply to fill in the residuals for the final term
2907  for(unsigned i=0;i<raw_ndof;i++)
2908  {
2909  unsigned local_eqn = elem_pt->eqn_number(i);
2910  //Add the slack parameter to the final residuals
2911  residuals[i] += Sigma*Psi.global_value(local_eqn)/
2912  Count.global_value(local_eqn);
2913  }
2914  }
2915  break;
2916 
2917  //If we are solving the augmented-by-one system
2918  case Block_augmented_J:
2919  {
2920  //Get the basic residuals
2921  elem_pt->get_residuals(residuals);
2922 
2923  //Zero the final residual
2924  residuals[raw_ndof] = 0.0;
2925  //Now multiply to fill in the residuals for the final term
2926  for(unsigned i=0;i<raw_ndof;i++)
2927  {
2928  unsigned local_eqn = elem_pt->eqn_number(i);
2929  //Add the slack parameter to the final residuals
2930  residuals[i] += Sigma*Psi.global_value(local_eqn)/
2931  Count.global_value(local_eqn);
2932  //Final term that specifies the symmetry
2933  residuals[raw_ndof] +=
2934  ((*Problem_pt->global_dof_pt(local_eqn))*Psi.global_value(local_eqn))/
2935  Count.global_value(local_eqn);
2936  }
2937  }
2938  break;
2939 
2940  //Otherwise we are solving the fully augemented system
2941  case Full_augmented:
2942  {
2943  DenseMatrix<double> jacobian(raw_ndof);
2944 
2945  //Get the basic residuals and jacobian
2946  elem_pt->get_jacobian(residuals,jacobian);
2947 
2948  //Initialise the final residuals
2949  residuals[raw_ndof] = 0.0;
2950  residuals[2*raw_ndof+1] = -1.0/Nelement;
2951 
2952  //Now multiply to fill in the residuals associated
2953  //with the null vector condition
2954  for(unsigned i=0;i<raw_ndof;i++)
2955  {
2956  unsigned local_eqn = elem_pt->eqn_number(i);
2957  residuals[raw_ndof+1+i] = 0.0;
2958  for(unsigned j=0;j<raw_ndof;j++)
2959  {
2960  unsigned local_unknown = elem_pt->eqn_number(j);
2961  residuals[raw_ndof+1+i] += jacobian(i,j)*Y.global_value(local_unknown);
2962  }
2963 
2964  //Add the slack parameter to the governing equations
2965  residuals[i] += Sigma*Psi.global_value(local_eqn)/
2966  Count.global_value(local_eqn);
2967 
2968  //Specify the symmetry
2969  residuals[raw_ndof] +=
2970  ((*Problem_pt->global_dof_pt(local_eqn))*Psi.global_value(local_eqn))/
2971  Count.global_value(local_eqn);
2972  //Specify the normalisation
2973  residuals[2*raw_ndof+1] += (Y.global_value(local_eqn)*
2974  C.global_value(local_eqn))/
2975  Count.global_value(local_eqn);
2976  }
2977  }
2978  break;
2979 
2980  default:
2981  std::ostringstream error_stream;
2982  error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
2983  << " not " << Solve_which_system << "\n";
2984  throw OomphLibError(error_stream.str(),
2985  OOMPH_CURRENT_FUNCTION,
2986  OOMPH_EXCEPTION_LOCATION);
2987  }
2988 
2989  }
2990 
2991  //======================================================================
2992  /// \short Calculate the elemental Jacobian matrix "d equation
2993  /// / d variable".
2994  //======================================================================
2996  Vector<double> &residuals,
2997  DenseMatrix<double> &jacobian)
2998  {
2999 
3000  unsigned augmented_ndof = ndof(elem_pt);
3001  unsigned raw_ndof = elem_pt->ndof();
3002 
3003  //Which system are we solving
3004  switch(Solve_which_system)
3005  {
3006  //If we are solving the original system
3007  case Block_J:
3008  {
3009  //get the raw residuals and jacobian
3010  elem_pt->get_jacobian(residuals,jacobian);
3011 
3012  //Now multiply to fill in the residuals for the final term
3013  for(unsigned i=0;i<raw_ndof;i++)
3014  {
3015  unsigned local_eqn = elem_pt->eqn_number(i);
3016  //Add the slack parameter to the final residuals
3017  residuals[i] += Sigma*Psi.global_value(local_eqn)/
3018  Count.global_value(local_eqn);
3019  }
3020  }
3021  break;
3022 
3023  //If we are solving the augmented-by-one system
3024  case Block_augmented_J:
3025  {
3026  //Get the full residuals, we need them
3027  get_residuals(elem_pt,residuals);
3028 
3029  //Need to get the raw jacobian (and raw residuals)
3030  Vector<double> newres(augmented_ndof);
3031  elem_pt->get_jacobian(newres,jacobian);
3032 
3033  //Now do finite differencing stuff
3034  const double FD_step = 1.0e-8;
3035  //Fill in the first lot of finite differences
3036  {
3037  //increase the global parameter
3038  double *unknown_pt = Parameter_pt;
3039  double init = *unknown_pt;
3040  *unknown_pt += FD_step;
3041 
3042  //Not do any possible updates
3043  Problem_pt
3045 
3046  //Get the new (modified) residuals
3047  get_residuals(elem_pt,newres);
3048 
3049  //The final column is given by the difference
3050  //between the residuals
3051  for(unsigned n=0;n<raw_ndof;n++)
3052  {
3053  jacobian(n,augmented_ndof-1) = (newres[n] - residuals[n])/FD_step;
3054  }
3055  //Reset the global parameter
3056  *unknown_pt = init;
3057 
3058  //Now do any possible updates
3059  Problem_pt
3061  }
3062 
3063  //Fill in the bottom row
3064  for(unsigned n=0;n<raw_ndof;n++)
3065  {
3066  unsigned local_eqn = elem_pt->eqn_number(n);
3067  jacobian(augmented_ndof-1,n)
3068  = Psi.global_value(local_eqn)/Count.global_value(local_eqn);
3069  }
3070  }
3071  break;
3072 
3073  //Otherwise solving the full system
3074  case Full_augmented:
3075  {
3076  /// ALICE:
3078  //Get the basic jacobian and residuals
3079  elem_pt->get_jacobian(residuals,jacobian);
3080  //get the proper residuals
3081  get_residuals(elem_pt,residuals);
3082 
3083  //Now fill in the next diagonal jacobian entry
3084  for(unsigned n=0;n<raw_ndof;n++)
3085  {
3086  for(unsigned m=0;m<raw_ndof;m++)
3087  {
3088  jacobian(raw_ndof+1+n,raw_ndof+1+m) = jacobian(n,m);
3089  }
3090  unsigned local_eqn = elem_pt->eqn_number(n);
3091  //Add in the sigma contribution
3092  jacobian(n,2*raw_ndof+1) = Psi.global_value(local_eqn)/
3093  Count.global_value(local_eqn);
3094  //Symmetry constraint
3095  jacobian(raw_ndof,n) = Psi.global_value(local_eqn)/
3096  Count.global_value(local_eqn);
3097  //Non-zero constraint
3098  jacobian(2*raw_ndof+1,raw_ndof+1+n) = C.global_value(local_eqn)/
3099  Count.global_value(local_eqn);
3100  }
3101 
3102  //Finite difference the remaining blocks
3103  const double FD_step = 1.0e-8;
3104 
3105  Vector<double> newres_p(augmented_ndof);
3106 
3107  //Loop over the ndofs only
3108  for(unsigned n=0;n<raw_ndof;++n)
3109  {
3110  //Get the original (unaugmented) global equation number
3111  unsigned long global_eqn = Assembly_handler_pt->eqn_number(elem_pt,n);
3112  double* unknown_pt = Problem_pt->global_dof_pt(global_eqn);
3113  double init = *unknown_pt;
3114  *unknown_pt += FD_step;
3116  //Get the new residuals
3117  get_residuals(elem_pt,newres_p);
3118  //Fill in the entries in the block d(Jy)/dx
3119  for(unsigned m=0;m<raw_ndof;m++)
3120  {
3121  jacobian(raw_ndof+1+m,n) =
3122  (newres_p[raw_ndof+1+m] - residuals[raw_ndof+1+m])/(FD_step);
3123  }
3124 
3125  //Reset the unknown
3126  *unknown_pt = init;
3128  }
3129 
3130  {
3131  //Now increase the global parameter
3132  double* unknown_pt = Parameter_pt;
3133  double init = *unknown_pt;
3134  *unknown_pt += FD_step;
3135 
3137  //Get the new residuals
3138  get_residuals(elem_pt,newres_p);
3139 
3140  //Add in the first block d/dx
3141  for(unsigned m=0;m<raw_ndof;m++)
3142  {
3143  jacobian(m,raw_ndof) =
3144  (newres_p[m] - residuals[m])/FD_step;
3145  }
3146  //Add in the second block d/dy
3147  for(unsigned m=raw_ndof+1;m<augmented_ndof-1;m++)
3148  {
3149  jacobian(m,raw_ndof) =
3150  (newres_p[m] - residuals[m])/FD_step;
3151  }
3152 
3153  //Reset the unknown
3154  *unknown_pt = init;
3156  }
3157  }
3158  break;
3159 
3160  default:
3161  std::ostringstream error_stream;
3162  error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
3163  << " not " << Solve_which_system << "\n";
3164  throw OomphLibError(error_stream.str(),
3165  OOMPH_CURRENT_FUNCTION,
3166  OOMPH_EXCEPTION_LOCATION);
3167  }
3168  /// ALICE
3169  Problem_pt
3171  }
3172 
3173 
3174 
3175  //==============================================================
3176  ///Get the derivatives of the residuals with respect to a parameter
3177  //==============================================================
3179  GeneralisedElement* const &elem_pt, double* const &parameter_pt,
3180  Vector<double> &dres_dparam)
3181  {
3182  //Need to get raw residuals and jacobian
3183  unsigned raw_ndof = elem_pt->ndof();
3185  //Find out which system we are solving
3186  switch(Solve_which_system)
3187  {
3188  //If we are solving the original system
3189  case Block_J:
3190  {
3191  //get the basic residual derivatives
3192  elem_pt->get_dresiduals_dparameter(parameter_pt,dres_dparam);
3193  //Slack parameter term does not depened explicitly on the parameter
3194  //so is not added
3195  }
3196  break;
3197 
3198  //If we are solving the augmented-by-one system
3199  case Block_augmented_J:
3200  {
3201  //Get the basic residual derivatives
3202  elem_pt->get_dresiduals_dparameter(parameter_pt,dres_dparam);
3203 
3204  //Zero the final residuals derivative
3205  dres_dparam[raw_ndof] = 0.0;
3206 
3207  //Other terms must not depend on the parameter
3208  }
3209  break;
3210 
3211  //Otherwise we are solving the fully augemented system
3212  case Full_augmented:
3213  {
3214  DenseMatrix<double> djac_dparam(raw_ndof);
3215 
3216  //Get the basic residuals and jacobian derivatives
3217  elem_pt->get_djacobian_dparameter(parameter_pt,dres_dparam,djac_dparam);
3218 
3219  //The "final" residual derivatives do not depend on the parameter
3220  dres_dparam[raw_ndof] = 0.0;
3221  dres_dparam[2*raw_ndof+1] = 0.0;
3222 
3223  //Now multiply to fill in the residuals associated
3224  //with the null vector condition
3225  for(unsigned i=0;i<raw_ndof;i++)
3226  {
3227  dres_dparam[raw_ndof+1+i] = 0.0;
3228  for(unsigned j=0;j<raw_ndof;j++)
3229  {
3230  unsigned local_unknown = elem_pt->eqn_number(j);
3231  dres_dparam[raw_ndof+1+i] +=
3232  djac_dparam(i,j)*Y.global_value(local_unknown);
3233  }
3234  }
3235  }
3236  break;
3237 
3238  default:
3239  std::ostringstream error_stream;
3240  error_stream << "The Solve_which_system flag can only take values 0, 1, 2"
3241  << " not " << Solve_which_system << "\n";
3242  throw OomphLibError(error_stream.str(),
3243  OOMPH_CURRENT_FUNCTION,
3244  OOMPH_EXCEPTION_LOCATION);
3245  }
3246 
3247  }
3248 
3249 
3250 
3251  //========================================================================
3252  /// Overload the derivative of the residuals and jacobian
3253  /// with respect to a parameter so that it breaks because it should not
3254  /// be required
3255  //========================================================================
3257  GeneralisedElement* const &elem_pt,
3258  double* const &parameter_pt,
3259  Vector<double> &dres_dparam,
3260  DenseMatrix<double> &djac_dparam)
3261  {
3262  std::ostringstream error_stream;
3263  error_stream <<
3264  "This function has not been implemented because it is not required\n";
3265  error_stream << "in standard problems.\n";
3266  error_stream <<
3267  "If you find that you need it, you will have to implement it!\n\n";
3268 
3269  throw OomphLibError(error_stream.str(),
3270  OOMPH_CURRENT_FUNCTION,
3271  OOMPH_EXCEPTION_LOCATION);
3272  }
3273 
3274 
3275  //=====================================================================
3276  /// Overload the hessian vector product function so that
3277  /// it calls the underlying element's hessian function.
3278  //========================================================================
3280  GeneralisedElement* const &elem_pt,
3281  Vector<double> const &Y,
3282  DenseMatrix<double> const &C,
3283  DenseMatrix<double> &product)
3284  {
3285  elem_pt->get_hessian_vector_products(Y,C,product);
3286  }
3287 
3288 
3289 
3290  //==========================================================================
3291  /// Return the eigenfunction(s) associated with the bifurcation that
3292  /// has been detected in bifurcation tracking problems
3293  //==========================================================================
3295  Vector<DoubleVector> &eigenfunction)
3296  {
3297  //There is only one (real) null vector
3298  eigenfunction.resize(1);
3299  //Rebuild the vector
3300  eigenfunction[0].build(this->Dof_distribution_pt,0.0);
3301  //Set the value to be the null vector
3302  const unsigned n_row_local = eigenfunction[0].nrow_local();
3303  for(unsigned n=0;n<n_row_local;n++)
3304  {
3305  eigenfunction[0][n] = Y[n];
3306  }
3307  }
3308 
3309 #ifdef OOMPH_HAS_MPI
3310 //=====================================================================
3311 //Synchronise the required data
3312 //====================================================================
3314  {
3315  //Only need to bother if the problem is distributed
3316  if(Distributed)
3317  {
3318  //Need only to synchronise the eigenfunction
3319  Y.synchronise();
3320  //Also need to synchronise the parameter and the slack parameter
3321  double broadcast_data[2];
3322  broadcast_data[0] = *Parameter_pt;
3323  broadcast_data[1] = Sigma;
3324  //Broadcast from the root to all processors
3325  MPI_Bcast(broadcast_data,2,MPI_DOUBLE,0,
3326  Problem_pt->communicator_pt()->mpi_comm());
3327 
3328  //Now copy received the values back
3329  *Parameter_pt = broadcast_data[0];
3330  Sigma = broadcast_data[1];
3331  }
3332  }
3333 #endif
3334 
3335  //====================================================================
3336  /// Set to solve the augmented-by-one block system.
3337  //===================================================================
3339  {
3340  //Only bother to do anything if we haven't already set the flag
3342  {
3343  //If we were solving the system with the original jacobian add the
3344  //parameter back
3346  {Problem_pt->Dof_pt.push_back(Parameter_pt);}
3347 
3348  //Restrict the problem to the standard variables and
3349  //the bifurcation parameter only
3351  (Problem_pt->communicator_pt(),Ndof+1,false);
3352  Problem_pt->Dof_pt.resize(Ndof+1);
3353 
3354 //Problem_pt->Dof_pt.resize(Ndof+1);
3355 // Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
3356 // Ndof+1,false);
3357 
3358  //Remove all previous sparse storage used during Jacobian assembly
3360 
3361  //Set the solve flag
3363  }
3364  }
3365 
3366 
3367 
3368  //==============================================================
3369  /// Set to solve the block system. The boolean flag is used
3370  /// to specify whether the block decomposition should use exactly
3371  /// the same jacobian as the original system.
3372  //==============================================================
3374  {
3375  //Only bother to do anything if we haven't already set the flag
3377  {
3378  //Restrict the problem to the standard variables
3379  Problem_pt->Dof_pt.resize(this->Dof_distribution_pt->nrow_local());
3381 
3382  //Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
3383  //Ndof,false);
3384 
3385  //Remove all previous sparse storage used during Jacobian assembly
3387 
3388  //Set the solve flag
3390  }
3391  }
3392 
3393  //==============================================================
3394  /// \short Solve non-block system
3395  //==============================================================
3397  {
3398  //Only do something if we are not solving the full system
3400  {
3401 #ifdef OOMPH_HAS_MPI
3402  int my_rank = Problem_pt->communicator_pt()->my_rank();
3403 #endif
3404  //If we were solving the problem without any augementation
3405  //add the parameter again
3407  {
3408 #ifdef OOMPH_HAS_MPI
3409 
3410  if((!Distributed) || (my_rank==0))
3411 #endif
3412  {
3413  Problem_pt->Dof_pt.push_back(Parameter_pt);
3414  }
3415  }
3416 
3417  //Always add the additional unknowns back into the problem
3418  const unsigned n_dof_local = Dof_distribution_pt->nrow_local();
3419  for(unsigned n=0;n<n_dof_local;n++)
3420  {
3421  Problem_pt->Dof_pt.push_back(&Y[n]);
3422  }
3423 
3424 
3425 
3426 #ifdef OOMPH_HAS_MPI
3427  if((!Distributed) || (my_rank==0))
3428 #endif
3429  //Add the slack parameter to the problem on the first processor
3430  //if distributed
3431  {
3432  Problem_pt->Dof_pt.push_back(&Sigma);
3433  }
3434 
3435 
3436  //Delete the distribtion created in the augmented block solve
3437  if(Problem_pt->Dof_distribution_pt!=this->Dof_distribution_pt)
3438  {delete Problem_pt->Dof_distribution_pt;}
3439 
3440  // update the Dof distribution
3443  //Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
3444  // Ndof*2+2,false);
3445  //Remove all previous sparse storage used during Jacobian assembly
3447 
3448  //Set the solve flag
3450  }
3451  }
3452 
3453 
3454  ///////////////////////////////////////////////////////////////////////
3455  // Non-inline functions for the BlockHopfLinearSolver class
3456  //////////////////////////////////////////////////////////////////////
3457 
3458 
3459  //======================================================================
3460  /// Clean up the memory that may have been allocated by the solver
3461  //=====================================================================
3463  {
3464  if(A_pt!=0) {delete A_pt;}
3465  if(E_pt!=0) {delete E_pt;}
3466  if(G_pt!=0) {delete G_pt;}
3467  }
3468 
3469  //===================================================================
3470  /// Use a block factorisation to solve the augmented system
3471  /// associated with a Hopf bifurcation.
3472  //===================================================================
3473  void BlockHopfLinearSolver::solve(Problem* const &problem_pt,
3474  DoubleVector &result)
3475  {
3476 
3477  // if the result is setup then it should not be distributed
3478 #ifdef PARANOID
3479  if (result.built())
3480  {
3481  if (result.distributed())
3482  {
3483  throw OomphLibError("The result vector must not be distributed",
3484  OOMPH_CURRENT_FUNCTION,
3485  OOMPH_EXCEPTION_LOCATION);
3486  }
3487  }
3488 #endif
3489 
3490  //Locally cache the pointer to the handler.
3491  HopfHandler* handler_pt =
3492  static_cast<HopfHandler*>(problem_pt->assembly_handler_pt());
3493 
3494  //Find the number of dofs in the augmented problem
3495  unsigned n_dof = problem_pt->ndof();
3496 
3497  // create the linear algebra distribution for this solver
3498  // currently only global (non-distributed) distributions are allowed
3499  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
3500  this->build_distribution(dist);
3501 
3502  //Firstly, let's calculate the derivative of the residuals wrt
3503  //the parameter
3504  DoubleVector dRdparam(this->distribution_pt(),0.0);
3505 
3506  const double FD_step = 1.0e-8;
3507  {
3508  //Cache pointer to the parameter (second last entry in the vector
3509  double *param_pt = &problem_pt->dof(n_dof-2);
3510  //Backup the parameter
3511  double old_var = *param_pt;;
3512  //Increment the parameter
3513  *param_pt += FD_step;
3515 
3516  //Now calculate the new residuals
3517  problem_pt->get_residuals(dRdparam);
3518 
3519  //Now calculate the difference assume original residuals in resul
3520  for(unsigned n=0;n<n_dof;n++)
3521  {
3522  dRdparam[n] = (dRdparam[n] - result[n])/FD_step;
3523  }
3524 
3525  //Reset the parameter
3526  *param_pt = old_var;
3528  }
3529 
3530  //Switch the handler to "standard" mode
3531  handler_pt->solve_standard_system();
3532 
3533  //Find out the number of dofs in the non-standard problem
3534  n_dof = problem_pt->ndof();
3535 
3536  // update the distribution pt
3537  dist.build(problem_pt->communicator_pt(),n_dof,false);
3538  this->build_distribution(dist);
3539 
3540  //Setup storage for temporary vector
3541  DoubleVector y1(this->distribution_pt(),0.0),
3542  alpha(this->distribution_pt(),0.0);
3543 
3544  //Allocate storage for A which can be used in the resolve
3545  if(A_pt!=0) {delete A_pt;}
3546  A_pt = new DoubleVector(this->distribution_pt(),0.0);
3547 
3548  //We are going to do resolves using the underlying linear solver
3550 
3551  //Solve the first system Jy1 = R
3552  Linear_solver_pt->solve(problem_pt,y1);
3553 
3554  //This should have set the sign of the jacobian, store it
3555  int sign_of_jacobian = problem_pt->sign_of_jacobian();
3556 
3557  //Now let's get the appropriate bit of alpha
3558  for(unsigned n=0;n<n_dof;n++) {alpha[n] = dRdparam[n];}
3559 
3560  //Resolve to find A
3561  Linear_solver_pt->resolve(alpha,*A_pt);
3562 
3563  //Now set to the complex system
3564  handler_pt->solve_complex_system();
3565 
3566  // update the distribution
3567  dist.build(problem_pt->communicator_pt(),n_dof*2,false);
3568  this->build_distribution(dist);
3569 
3570  //Resize the stored vector G
3571  if(G_pt!=0) {delete G_pt;}
3572  G_pt = new DoubleVector(this->distribution_pt(),0.0);
3573 
3574  //Solve the first Zg = (Mz -My)
3575  Linear_solver_pt->solve(problem_pt,*G_pt);
3576 
3577  //This should have set the sign of the complex matrix's determinant, multiply
3578  sign_of_jacobian *= problem_pt->sign_of_jacobian();
3579 
3580  //We can now construct our multipliers
3581  //Prepare to scale
3582  double dof_length=0.0, a_length=0.0, y1_length=0.0;
3583  //Loop over the standard number of dofs
3584  for(unsigned n=0;n<n_dof;n++)
3585  {
3586  if(std::fabs(problem_pt->dof(n)) > dof_length)
3587  {dof_length = std::fabs(problem_pt->dof(n));}
3588  if(std::fabs((*A_pt)[n]) > a_length) {a_length = std::fabs((*A_pt)[n]);}
3589  if(std::fabs(y1[n]) > y1_length)
3590  {y1_length = std::fabs(y1[n]);}
3591  }
3592 
3593  double a_mult = dof_length/a_length;
3594  double y1_mult = dof_length/y1_length;
3595  a_mult += FD_step; y1_mult += FD_step;
3596  a_mult *= FD_step; y1_mult *= FD_step;
3597 
3598  //Local storage for the product terms
3599  Vector<double> Jprod_a(2*n_dof,0.0), Jprod_y1(2*n_dof,0.0);
3600 
3601  //Temporary storage
3602  Vector<double> rhs(2*n_dof);
3603 
3604  //find the number of elements
3605  unsigned long n_element = problem_pt->mesh_pt()->nelement();
3606 
3607  //Calculate the product of the jacobian matrices, etc
3608  for(unsigned long e=0;e<n_element;e++)
3609  {
3610  GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
3611  //Loop over the ndofs in each element
3612  unsigned n_var = elem_pt->ndof();
3613  //Get the jacobian matrices
3614  DenseMatrix<double> jac(n_var), jac_a(n_var), jac_y1(n_var);
3615  DenseMatrix<double> M(n_var), M_a(n_var), M_y1(n_var);
3616  //Get unperturbed jacobian
3617  elem_pt->get_jacobian_and_mass_matrix(rhs,jac,M);
3618 
3619  //Backup the dofs
3620  Vector<double> dof_bac(n_var);
3621  //Perturb the original dofs
3622  for(unsigned n=0;n<n_var;n++)
3623  {
3624  unsigned eqn_number = elem_pt->eqn_number(n);
3625  dof_bac[n] = problem_pt->dof(eqn_number);
3626  //Pertub by vector A
3627  problem_pt->dof(eqn_number) += a_mult*(*A_pt)[eqn_number];
3628  }
3629 
3630  //Now get the new jacobian and mass matrix
3631  elem_pt->get_jacobian_and_mass_matrix(rhs,jac_a,M_a);
3632 
3633  //Perturb the dofs
3634  for(unsigned n=0;n<n_var;n++)
3635  {
3636  unsigned eqn_number = elem_pt->eqn_number(n);
3637  problem_pt->dof(eqn_number) = dof_bac[n];
3638  //Pertub by vector y1
3639  problem_pt->dof(eqn_number) += y1_mult*y1[eqn_number];
3640  }
3641 
3642  //Now get the new jacobian and mass matrix
3643  elem_pt->get_jacobian_and_mass_matrix(rhs,jac_y1,M_y1);
3644 
3645  //Reset the dofs
3646  for(unsigned n=0;n<n_var;n++)
3647  {
3648  unsigned eqn_number = elem_pt->eqn_number(n);
3649  problem_pt->dof(eqn_number) = dof_bac[n];
3650  }
3651 
3652  //OK, now work out the products
3653  for(unsigned n=0;n<n_var;n++)
3654  {
3655  unsigned eqn_number = elem_pt->eqn_number(n);
3656  double prod_a1=0.0, prod_y11=0.0;
3657  double prod_a2=0.0, prod_y12=0.0;
3658  for(unsigned m=0;m<n_var;m++)
3659  {
3660  const unsigned unknown = elem_pt->eqn_number(m);
3661  const double y = handler_pt->Phi[unknown];
3662  const double z = handler_pt->Psi[unknown];
3663  const double omega = handler_pt->Omega;
3664  //Real part (first line)
3665  prod_a1 += (jac_a(n,m) - jac(n,m))*y
3666  + omega*(M_a(n,m) - M(n,m))*z;
3667  prod_y11 += (jac_y1(n,m) - jac(n,m))*y
3668  + omega*(M_y1(n,m) - M(n,m))*z;
3669  //Imag par (second line)
3670  prod_a2 += (jac_a(n,m) - jac(n,m))*z
3671  - omega*(M_a(n,m) - M(n,m))*y;
3672  prod_y12 += (jac_y1(n,m) - jac(n,m))*z
3673  - omega*(M_y1(n,m) - M(n,m))*y;
3674  }
3675  Jprod_a[eqn_number] += prod_a1/a_mult;
3676  Jprod_y1[eqn_number] += prod_y11/y1_mult;
3677  Jprod_a[n_dof + eqn_number] += prod_a2/a_mult;
3678  Jprod_y1[n_dof + eqn_number] += prod_y12/y1_mult;
3679  }
3680  }
3681 
3682  //The assumption here is still that the result has been set to the
3683  //residuals.
3684  for(unsigned n=0;n<2*n_dof;n++)
3685  {
3686  rhs[n] = result[n_dof+n] - Jprod_y1[n];
3687  }
3688 
3689  //Temporary storage
3690  DoubleVector y2(this->distribution_pt(),0.0);
3691 
3692  // DoubleVector for rhs
3693  DoubleVector rhs2(this->distribution_pt(),0.0);
3694  for (unsigned i = 0; i < 2*n_dof; i++)
3695  {
3696  rhs2[i] = rhs[i];
3697  }
3698 
3699  //Solve it
3700  Linear_solver_pt->resolve(rhs2,y2);
3701 
3702  //Assemble the next RHS
3703  for(unsigned n=0;n<2*n_dof;n++)
3704  {
3705  rhs[n] = dRdparam[n_dof+n] - Jprod_a[n];
3706  }
3707 
3708  //Resive the storage
3709  if(E_pt!=0) {delete E_pt;}
3710  E_pt = new DoubleVector(this->distribution_pt(),0.0);
3711 
3712  //Solve for the next RHS
3713  for (unsigned i = 0; i < 2*n_dof; i++)
3714  {
3715  rhs2[i] = rhs[i];
3716  }
3717  Linear_solver_pt->resolve(rhs2,*E_pt);
3718 
3719  //We can now calculate the final corrections
3720  //We need to work out a large number of dot products
3721  double dot_c = 0.0, dot_d = 0.0, dot_e = 0.0, dot_f = 0.0, dot_g=0.0;
3722  double dot_h = 0.0;
3723 
3724  for(unsigned n=0;n<n_dof;n++)
3725  {
3726  //Get the appopriate entry
3727  const double Cn = handler_pt->C[n];
3728  dot_c += Cn*y2[n];
3729  dot_d += Cn*y2[n_dof+n];
3730  dot_e += Cn*(*E_pt)[n];
3731  dot_f += Cn*(*E_pt)[n_dof+n];
3732  dot_g += Cn*(*G_pt)[n];
3733  dot_h += Cn*(*G_pt)[n_dof+n];
3734  }
3735 
3736  //Now we should be able to work out the corrections
3737  double denom = dot_e*dot_h - dot_g*dot_f;
3738 
3739  //Copy the previous residuals
3740  double R31 = result[3*n_dof], R32 = result[3*n_dof+1];
3741  //Delta parameter
3742  const double delta_param =
3743  ((R32 - dot_d)*dot_g - (R31 - dot_c)*dot_h)/denom;
3744  //Delta frequency
3745  const double delta_w = -((R32 - dot_d) + dot_f*delta_param)/(dot_h);
3746 
3747  //Load into the result vector
3748  result[3*n_dof] = delta_param;
3749  result[3*n_dof+1] = delta_w;
3750 
3751  //The corrections to the null vector
3752  for(unsigned n=0;n<2*n_dof;n++)
3753  {
3754  result[n_dof + n] = y2[n] - (*E_pt)[n]*delta_param - (*G_pt)[n]*delta_w;
3755  }
3756 
3757  //Finally add the corrections to the unknowns
3758  for(unsigned n=0;n<n_dof;n++)
3759  {
3760  result[n] = y1[n] - (*A_pt)[n]*delta_param;
3761  }
3762 
3763  //The sign of the jacobian is the previous signs multiplied by the
3764  //sign of the denominator
3765  problem_pt->sign_of_jacobian() = sign_of_jacobian*
3766  static_cast<int>(std::fabs(denom)/denom);
3767 
3768  //Switch things to our full solver
3769  handler_pt->solve_full_system();
3770 
3771  //If we are not storing things, clear the results
3772  if(!Enable_resolve)
3773  {
3774  //We no longer need to store the matrix
3776  delete A_pt; A_pt=0;
3777  delete E_pt; E_pt=0;
3778  delete G_pt; G_pt=0;
3779  }
3780  //Otherwise, also store the problem pointer
3781  else
3782  {
3783  Problem_pt = problem_pt;
3784  }
3785  }
3786 
3787  //=====================================================================
3788  /// Solve for two right hand sides, required for (efficient) continuation
3789  /// because otherwise we have to store the inverse of the jacobian
3790  /// and the complex equivalent ... nasty.
3791  //=====================================================================
3793  DoubleVector &result,
3794  const DoubleVector &rhs2,
3795  DoubleVector &result2)
3796  {
3797 
3798  // if the result is setup then it should not be distributed
3799 #ifdef PARANOID
3800  if (result.built())
3801  {
3802  if (result.distributed())
3803  {
3804  throw OomphLibError("The result vector must not be distributed",
3805  OOMPH_CURRENT_FUNCTION,
3806  OOMPH_EXCEPTION_LOCATION);
3807  }
3808  }
3809  if (result2.built())
3810  {
3811  if (result2.distributed())
3812  {
3813  throw OomphLibError("The result2 vector must not be distributed",
3814  OOMPH_CURRENT_FUNCTION,
3815  OOMPH_EXCEPTION_LOCATION);
3816  }
3817  }
3818 #endif
3819 
3820  //Locally cache the pointer to the handler.
3821  HopfHandler* handler_pt =
3822  static_cast<HopfHandler*>(problem_pt->assembly_handler_pt());
3823 
3824  //Find the number of dofs in the augmented problem
3825  unsigned n_dof = problem_pt->ndof();
3826 
3827  // create the linear algebra distribution for this solver
3828  // currently only global (non-distributed) distributions are allowed
3829  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
3830  this->build_distribution(dist);
3831 
3832  // if the result vector is not setup then rebuild with distribution = global
3833  if (!result.built())
3834  {
3835  result.build(this->distribution_pt(),0.0);
3836  }
3837  if (!result2.built())
3838  {
3839  result2.build(this->distribution_pt(),0.0);
3840  }
3841 
3842  //Firstly, let's calculate the derivative of the residuals wrt
3843  //the parameter
3844  DoubleVector dRdparam(this->distribution_pt(),0.0);
3845 
3846  const double FD_step = 1.0e-8;
3847  {
3848  //Cache pointer to the parameter (second last entry in the vector
3849  double *param_pt = &problem_pt->dof(n_dof-2);
3850  //Backup the parameter
3851  double old_var = *param_pt;;
3852  //Increment the parameter
3853  *param_pt += FD_step;
3855 
3856  //Now calculate the new residuals
3857  problem_pt->get_residuals(dRdparam);
3858 
3859  //Now calculate the difference assume original residuals in resul
3860  for(unsigned n=0;n<n_dof;n++)
3861  {
3862  dRdparam[n] = (dRdparam[n] - result[n])/FD_step;
3863  }
3864 
3865  //Reset the parameter
3866  *param_pt = old_var;
3868  }
3869 
3870  //Switch the handler to "standard" mode
3871  handler_pt->solve_standard_system();
3872 
3873  //Find out the number of dofs in the non-standard problem
3874  n_dof = problem_pt->ndof();
3875 
3876  //
3877  dist.build(problem_pt->communicator_pt(),n_dof,false);
3878  this->build_distribution(dist);
3879 
3880  //Setup storage for temporary vector
3881  DoubleVector y1(this->distribution_pt(),0.0),
3882  alpha(this->distribution_pt(),0.0),
3883  y1_resolve(this->distribution_pt(),0.0);
3884 
3885  //Allocate storage for A which can be used in the resolve
3886  if(A_pt!=0) {delete A_pt;}
3887  A_pt = new DoubleVector(this->distribution_pt(),0.0);
3888 
3889  //We are going to do resolves using the underlying linear solver
3891 
3892  //Solve the first system Jy1 =
3893  Linear_solver_pt->solve(problem_pt,y1);
3894 
3895  //This should have set the sign of the jacobian, store it
3896  int sign_of_jacobian = problem_pt->sign_of_jacobian();
3897 
3898  //Now let's get the appropriate bit of alpha
3899  for(unsigned n=0;n<n_dof;n++) {alpha[n] = dRdparam[n];}
3900 
3901  //Resolve to find A
3902  Linear_solver_pt->resolve(alpha,*A_pt);
3903 
3904  //Get the solution for the second rhs
3905  for(unsigned n=0;n<n_dof;n++) {alpha[n] = rhs2[n];}
3906 
3907  //Resolve to find y1_resolve
3908  Linear_solver_pt->resolve(alpha,y1_resolve);
3909 
3910  //Now set to the complex system
3911  handler_pt->solve_complex_system();
3912 
3913  // rebuild the Distribution
3914  dist.build(problem_pt->communicator_pt(),n_dof*2,false);
3915  this->build_distribution(dist);
3916 
3917  //Resize the stored vector G
3918  if(G_pt!=0) {delete G_pt;}
3919  G_pt = new DoubleVector(this->distribution_pt(),0.0);
3920 
3921  //Solve the first Zg = (Mz -My)
3922  Linear_solver_pt->solve(problem_pt,*G_pt);
3923 
3924  //This should have set the sign of the complex matrix's determinant, multiply
3925  sign_of_jacobian *= problem_pt->sign_of_jacobian();
3926 
3927  //We can now construct our multipliers
3928  //Prepare to scale
3929  double dof_length=0.0, a_length=0.0, y1_length=0.0, y1_resolve_length=0.0;
3930  //Loop over the standard number of dofs
3931  for(unsigned n=0;n<n_dof;n++)
3932  {
3933  if(std::fabs(problem_pt->dof(n)) > dof_length)
3934  {dof_length = std::fabs(problem_pt->dof(n));}
3935  if(std::fabs((*A_pt)[n]) > a_length) {a_length = std::fabs((*A_pt)[n]);}
3936  if(std::fabs(y1[n]) > y1_length)
3937  {y1_length = std::fabs(y1[n]);}
3938  if(std::fabs(y1_resolve[n]) > y1_resolve_length)
3939  {y1_resolve_length = std::fabs(y1[n]);}
3940  }
3941 
3942 
3943  double a_mult = dof_length/a_length;
3944  double y1_mult = dof_length/y1_length;
3945  double y1_resolve_mult = dof_length/y1_resolve_length;
3946  a_mult += FD_step; y1_mult += FD_step; y1_resolve_mult += FD_step;
3947  a_mult *= FD_step; y1_mult *= FD_step; y1_resolve_mult *= FD_step;
3948 
3949  //Local storage for the product terms
3950  Vector<double> Jprod_a(2*n_dof,0.0), Jprod_y1(2*n_dof,0.0);
3951  Vector<double> Jprod_y1_resolve(2*n_dof,0.0);
3952 
3953  //Temporary storage
3954  Vector<double> rhs(2*n_dof);
3955 
3956  //find the number of elements
3957  unsigned long n_element = problem_pt->mesh_pt()->nelement();
3958 
3959  //Calculate the product of the jacobian matrices, etc
3960  for(unsigned long e = 0;e<n_element;e++)
3961  {
3962  GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
3963  //Loop over the ndofs in each element
3964  unsigned n_var = elem_pt->ndof();
3965  //Get the jacobian matrices
3966  DenseMatrix<double> jac(n_var), jac_a(n_var), jac_y1(n_var),
3967  jac_y1_resolve(n_var);
3968  DenseMatrix<double> M(n_var), M_a(n_var), M_y1(n_var), M_y1_resolve(n_var);
3969  //Get unperturbed jacobian
3970  elem_pt->get_jacobian_and_mass_matrix(rhs,jac,M);
3971 
3972  //Backup the dofs
3973  Vector<double> dof_bac(n_var);
3974  //Perturb the original dofs
3975  for(unsigned n=0;n<n_var;n++)
3976  {
3977  unsigned eqn_number = elem_pt->eqn_number(n);
3978  dof_bac[n] = problem_pt->dof(eqn_number);
3979  //Pertub by vector A
3980  problem_pt->dof(eqn_number) += a_mult*(*A_pt)[eqn_number];
3981  }
3982 
3983  //Now get the new jacobian and mass matrix
3984  elem_pt->get_jacobian_and_mass_matrix(rhs,jac_a,M_a);
3985 
3986  //Perturb the dofs
3987  for(unsigned n=0;n<n_var;n++)
3988  {
3989  unsigned eqn_number = elem_pt->eqn_number(n);
3990  problem_pt->dof(eqn_number) = dof_bac[n];
3991  //Pertub by vector y1
3992  problem_pt->dof(eqn_number) += y1_mult*y1[eqn_number];
3993  }
3994 
3995  //Now get the new jacobian and mass matrix
3996  elem_pt->get_jacobian_and_mass_matrix(rhs,jac_y1,M_y1);
3997 
3998  //Perturb the dofs
3999  for(unsigned n=0;n<n_var;n++)
4000  {
4001  unsigned eqn_number = elem_pt->eqn_number(n);
4002  problem_pt->dof(eqn_number) = dof_bac[n];
4003  //Pertub by vector y1
4004  problem_pt->dof(eqn_number) += y1_resolve_mult*y1_resolve[eqn_number];
4005  }
4006 
4007  //Now get the new jacobian and mass matrix
4008  elem_pt->get_jacobian_and_mass_matrix(rhs,jac_y1_resolve,M_y1_resolve);
4009 
4010  //Reset the dofs
4011  for(unsigned n=0;n<n_var;n++)
4012  {
4013  unsigned eqn_number = elem_pt->eqn_number(n);
4014  problem_pt->dof(eqn_number) = dof_bac[n];
4015  }
4016 
4017  //OK, now work out the products
4018  for(unsigned n=0;n<n_var;n++)
4019  {
4020  unsigned eqn_number = elem_pt->eqn_number(n);
4021  double prod_a1=0.0, prod_y11=0.0, prod_y1_resolve1 = 0.0;
4022  double prod_a2=0.0, prod_y12=0.0, prod_y1_resolve2 = 0.0;
4023  for(unsigned m=0;m<n_var;m++)
4024  {
4025  const unsigned unknown = elem_pt->eqn_number(m);
4026  const double y = handler_pt->Phi[unknown];
4027  const double z = handler_pt->Psi[unknown];
4028  const double omega = handler_pt->Omega;
4029  //Real part (first line)
4030  prod_a1 += (jac_a(n,m) - jac(n,m))*y
4031  + omega*(M_a(n,m) - M(n,m))*z;
4032  prod_y11 += (jac_y1(n,m) - jac(n,m))*y
4033  + omega*(M_y1(n,m) - M(n,m))*z;
4034  prod_y1_resolve1 += (jac_y1_resolve(n,m) - jac(n,m))*y
4035  + omega*(M_y1_resolve(n,m) - M(n,m))*z;
4036  //Imag par (second line)
4037  prod_a2 += (jac_a(n,m) - jac(n,m))*z
4038  - omega*(M_a(n,m) - M(n,m))*y;
4039  prod_y12 += (jac_y1(n,m) - jac(n,m))*z
4040  - omega*(M_y1(n,m) - M(n,m))*y;
4041  prod_y1_resolve2 += (jac_y1_resolve(n,m) - jac(n,m))*z
4042  - omega*(M_y1_resolve(n,m) - M(n,m))*y;
4043  }
4044  Jprod_a[eqn_number] += prod_a1/a_mult;
4045  Jprod_y1[eqn_number] += prod_y11/y1_mult;
4046  Jprod_y1_resolve[eqn_number] += prod_y1_resolve1/y1_resolve_mult;
4047  Jprod_a[n_dof + eqn_number] += prod_a2/a_mult;
4048  Jprod_y1[n_dof + eqn_number] += prod_y12/y1_mult;
4049  Jprod_y1_resolve[n_dof + eqn_number] += prod_y1_resolve2/y1_resolve_mult;
4050  }
4051  }
4052 
4053  //The assumption here is still that the result has been set to the
4054  //residuals.
4055  for(unsigned n=0;n<2*n_dof;n++)
4056  {
4057  rhs[n] = result[n_dof+n] - Jprod_y1[n];
4058  }
4059 
4060  //Temporary storage
4061  DoubleVector y2(this->distribution_pt(),0.0);
4062 
4063  //Solve it
4064  DoubleVector temp_rhs(this->distribution_pt(),0.0);
4065  for (unsigned i = 0; i < 2*n_dof; i++)
4066  {
4067  temp_rhs[i] = rhs[i];
4068  }
4069  Linear_solver_pt->resolve(temp_rhs,y2);
4070 
4071  //Assemble the next RHS
4072  for(unsigned n=0;n<2*n_dof;n++)
4073  {
4074  rhs[n] = dRdparam[n_dof+n] - Jprod_a[n];
4075  }
4076 
4077  //Resive the storage
4078  if(E_pt!=0) {delete E_pt;}
4079  E_pt = new DoubleVector(this->distribution_pt(),0.0);
4080 
4081  //Solve for the next RHS
4082  for (unsigned i = 0; i < 2* n_dof; i++)
4083  {
4084  temp_rhs[i] = rhs[i];
4085  }
4086  Linear_solver_pt->resolve(temp_rhs,*E_pt);
4087 
4088  //Assemble the next RHS
4089  for(unsigned n=0;n<2*n_dof;n++)
4090  {
4091  rhs[n] = rhs2[n_dof+n] - Jprod_y1_resolve[n];
4092  }
4093 
4094  DoubleVector y2_resolve(this->distribution_pt(),0.0);
4095  for (unsigned i = 0; i < 2* n_dof; i++)
4096  {
4097  temp_rhs[i] = rhs[i];
4098  }
4099  Linear_solver_pt->resolve(temp_rhs,y2_resolve);
4100 
4101 
4102  //We can now calculate the final corrections
4103  //We need to work out a large number of dot products
4104  double dot_c = 0.0, dot_d = 0.0, dot_e = 0.0, dot_f = 0.0, dot_g=0.0;
4105  double dot_h = 0.0;
4106 
4107  double dot_c_resolve = 0.0, dot_d_resolve = 0.0;
4108 
4109  for(unsigned n=0;n<n_dof;n++)
4110  {
4111  //Get the appopriate entry
4112  const double Cn = handler_pt->C[n];
4113  dot_c += Cn*y2[n];
4114  dot_d += Cn*y2[n_dof+n];
4115  dot_c_resolve += Cn*y2_resolve[n];
4116  dot_d_resolve += Cn*y2_resolve[n_dof+n];
4117  dot_e += Cn*(*E_pt)[n];
4118  dot_f += Cn*(*E_pt)[n_dof+n];
4119  dot_g += Cn*(*G_pt)[n];
4120  dot_h += Cn*(*G_pt)[n_dof+n];
4121  }
4122 
4123  //Now we should be able to work out the corrections
4124  double denom = dot_e*dot_h - dot_g*dot_f;
4125 
4126  //Copy the previous residuals
4127  double R31 = result[3*n_dof], R32 = result[3*n_dof+1];
4128  //Delta parameter
4129  const double delta_param =
4130  ((R32 - dot_d)*dot_g - (R31 - dot_c)*dot_h)/denom;
4131  //Delta frequency
4132  const double delta_w = -((R32 - dot_d) + dot_f*delta_param)/(dot_h);
4133 
4134  //Corrections
4135  double R31_resolve = rhs2[3*n_dof], R32_resolve = rhs2[3*n_dof+1];
4136  //Delta parameter
4137  const double delta_param_resolve =
4138  ((R32_resolve - dot_d_resolve)*dot_g -
4139  (R31_resolve - dot_c_resolve)*dot_h)/denom;
4140  //Delta frequency
4141  const double delta_w_resolve =
4142  -((R32_resolve - dot_d_resolve) + dot_f*delta_param_resolve)/(dot_h);
4143 
4144 
4145  //Load into the result vector
4146  result[3*n_dof] = delta_param;
4147  result[3*n_dof+1] = delta_w;
4148 
4149  //The corrections to the null vector
4150  for(unsigned n=0;n<2*n_dof;n++)
4151  {
4152  result[n_dof + n] = y2[n] - (*E_pt)[n]*delta_param - (*G_pt)[n]*delta_w;
4153  }
4154 
4155  //Finally add the corrections to the unknowns
4156  for(unsigned n=0;n<n_dof;n++)
4157  {
4158  result[n] = y1[n] - (*A_pt)[n]*delta_param;
4159  }
4160 
4161  //Load into the result vector
4162  result2[3*n_dof] = delta_param_resolve;
4163  result2[3*n_dof+1] = delta_w_resolve;
4164 
4165  //The corrections to the null vector
4166  for(unsigned n=0;n<2*n_dof;n++)
4167  {
4168  result2[n_dof + n] = y2_resolve[n] - (*E_pt)[n]*delta_param_resolve
4169  - (*G_pt)[n]*delta_w_resolve;
4170  }
4171 
4172  //Finally add the corrections to the unknowns
4173  for(unsigned n=0;n<n_dof;n++)
4174  {
4175  result2[n] = y1_resolve[n] - (*A_pt)[n]*delta_param_resolve;
4176  }
4177 
4178 
4179  //The sign of the jacobian is the previous signs multiplied by the
4180  //sign of the denominator
4181  problem_pt->sign_of_jacobian() = sign_of_jacobian*
4182  static_cast<int>(std::fabs(denom)/denom);
4183 
4184  //Switch things to our full solver
4185  handler_pt->solve_full_system();
4186 
4187  //If we are not storing things, clear the results
4188  if(!Enable_resolve)
4189  {
4190  //We no longer need to store the matrix
4192  delete A_pt; A_pt=0;
4193  delete E_pt; E_pt=0;
4194  delete G_pt; G_pt=0;
4195  }
4196  //Otherwise, also store the problem pointer
4197  else
4198  {
4199  Problem_pt = problem_pt;
4200  }
4201 
4202  }
4203 
4204 
4205  //======================================================================
4206  //Hack the re-solve to use the block factorisation
4207  //======================================================================
4209  DoubleVector &result)
4210  {
4211  throw OomphLibError("resolve() is not implemented for this solver",
4212  OOMPH_CURRENT_FUNCTION,
4213  OOMPH_EXCEPTION_LOCATION);
4214 
4215  }
4216 
4217 
4218  ///////////////////////////////////////////////////////////////////////
4219  // Non-inline functions for the HopfHandler class
4220  //////////////////////////////////////////////////////////////////////
4221 
4222  //====================================================================
4223  /// Constructor: Initialise the hopf handler,
4224  /// by setting initial guesses for Phi, Psi and calculating Count.
4225  /// If the system changes, a new handler must be constructed.
4226  //===================================================================
4227  HopfHandler::HopfHandler(Problem* const &problem_pt,
4228  double* const &parameter_pt) :
4229  Solve_which_system(0), Parameter_pt(parameter_pt), Omega(0.0)
4230  {
4231  //Set the problem pointer
4232  Problem_pt = problem_pt;
4233  //Set the number of non-augmented degrees of freedom
4234  Ndof = problem_pt->ndof();
4235 
4236  // create the linear algebra distribution for this solver
4237  // currently only global (non-distributed) distributions are allowed
4238  LinearAlgebraDistribution* dist_pt = new
4239  LinearAlgebraDistribution(problem_pt->communicator_pt(),Ndof,false);
4240 
4241  //Resize the vectors of additional dofs
4242  Phi.resize(Ndof);
4243  Psi.resize(Ndof);
4244  C.resize(Ndof);
4245  Count.resize(Ndof,0);
4246 
4247  //Loop over all the elements in the problem
4248  unsigned n_element = problem_pt->mesh_pt()->nelement();
4249  for(unsigned e=0;e<n_element;e++)
4250  {
4251  GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
4252  //Loop over the local freedoms in an element
4253  unsigned n_var = elem_pt->ndof();
4254  for(unsigned n=0;n<n_var;n++)
4255  {
4256  //Increase the associated global equation number counter
4257  ++Count[elem_pt->eqn_number(n)];
4258  }
4259  }
4260 
4261  //Calculate the value Phi by
4262  //solving the system JPhi = dF/dlambda
4263 
4264  //Locally cache the linear solver
4265  LinearSolver* const linear_solver_pt = problem_pt->linear_solver_pt();
4266 
4267  //Save the status before entry to this routine
4268  bool enable_resolve = linear_solver_pt->is_resolve_enabled();
4269 
4270  //We need to do a resolve
4271  linear_solver_pt->enable_resolve();
4272 
4273  //Storage for the solution
4274  DoubleVector x(dist_pt,0.0);
4275 
4276  //Solve the standard problem, we only want to make sure that
4277  //we factorise the matrix, if it has not been factorised. We shall
4278  //ignore the return value of x.
4279  linear_solver_pt->solve(problem_pt,x);
4280 
4281  //Get the vector dresiduals/dparameter
4282  problem_pt->get_derivative_wrt_global_parameter(parameter_pt,x);
4283 
4284  // Copy rhs vector into local storage so it doesn't get overwritten
4285  // if the linear solver decides to initialise the solution vector, say,
4286  // which it's quite entitled to do!
4287  DoubleVector input_x(x);
4288 
4289  //Now resolve the system with the new RHS and overwrite the solution
4290  linear_solver_pt->resolve(input_x,x);
4291 
4292  //Restore the storage status of the linear solver
4293  if (enable_resolve)
4294  {
4295  linear_solver_pt->enable_resolve();
4296  }
4297  else
4298  {
4299  linear_solver_pt->disable_resolve();
4300  }
4301 
4302  //Normalise the solution x
4303  double length = 0.0;
4304  for(unsigned n=0;n<Ndof;n++) {length += x[n]*x[n];}
4305  length = sqrt(length);
4306 
4307  //Now add the real part of the null space components to the problem
4308  //unknowns and initialise it all
4309  //This is dumb at the moment ... fix with eigensolver?
4310  for(unsigned n=0;n<Ndof;n++)
4311  {
4312  problem_pt->Dof_pt.push_back(&Phi[n]);
4313  C[n] = Phi[n] = -x[n]/length;
4314  }
4315 
4316  //Set the imaginary part so that the appropriate residual is
4317  //zero initially (eigensolvers)
4318  for(unsigned n=0;n<Ndof;n+=2)
4319  {
4320  //Make sure that we are not at the end of an array of odd length
4321  if(n!=Ndof-1)
4322  {
4323  Psi[n] = C[n+1];
4324  Psi[n+1] = -C[n];
4325  }
4326  //If it's odd set the final entry to zero
4327  else
4328  {
4329  Psi[n] = 0.0;
4330  }
4331  }
4332 
4333  //Next add the imaginary parts of the null space components to the problem
4334  for(unsigned n=0;n<Ndof;n++)
4335  {
4336  problem_pt->Dof_pt.push_back(&Psi[n]);
4337  }
4338  //Now add the parameter
4339  problem_pt->Dof_pt.push_back(parameter_pt);
4340  //Finally add the frequency
4341  problem_pt->Dof_pt.push_back(&Omega);
4342 
4343  // rebuild the dof dist
4345  Ndof*3+2,false);
4346  //Remove all previous sparse storage used during Jacobian assembly
4348 
4349  // delete the dist_pt
4350  delete dist_pt;
4351  }
4352 
4353  //====================================================================
4354  /// Constructor: Initialise the hopf handler,
4355  /// by setting initial guesses for Phi, Psi, Omega and calculating Count.
4356  /// If the system changes, a new handler must be constructed.
4357  //===================================================================
4358  HopfHandler::HopfHandler(Problem* const &problem_pt,
4359  double* const &parameter_pt,
4360  const double &omega,
4361  const DoubleVector &phi,
4362  const DoubleVector &psi) :
4363  Solve_which_system(0), Parameter_pt(parameter_pt), Omega(omega)
4364  {
4365  //Set the problem pointer
4366  Problem_pt = problem_pt;
4367  //Set the number of non-augmented degrees of freedom
4368  Ndof = problem_pt->ndof();
4369 
4370  //Resize the vectors of additional dofs
4371  Phi.resize(Ndof);
4372  Psi.resize(Ndof);
4373  C.resize(Ndof);
4374  Count.resize(Ndof,0);
4375 
4376  //Loop over all the elements in the problem
4377  unsigned n_element = problem_pt->mesh_pt()->nelement();
4378  for(unsigned e=0;e<n_element;e++)
4379  {
4380  GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
4381  //Loop over the local freedoms in an element
4382  unsigned n_var = elem_pt->ndof();
4383  for(unsigned n=0;n<n_var;n++)
4384  {
4385  //Increase the associated global equation number counter
4386  ++Count[elem_pt->eqn_number(n)];
4387  }
4388  }
4389 
4390  //Normalise the guess for phi
4391  double length = 0.0;
4392  for(unsigned n=0;n<Ndof;n++) {length += phi[n]*phi[n];}
4393  length = sqrt(length);
4394 
4395  //Now add the real part of the null space components to the problem
4396  //unknowns and initialise it all
4397  for(unsigned n=0;n<Ndof;n++)
4398  {
4399  problem_pt->Dof_pt.push_back(&Phi[n]);
4400  C[n] = Phi[n] = phi[n]/length;
4401  Psi[n] = psi[n]/length;
4402  }
4403 
4404  //Next add the imaginary parts of the null space components to the problem
4405  for(unsigned n=0;n<Ndof;n++)
4406  {
4407  problem_pt->Dof_pt.push_back(&Psi[n]);
4408  }
4409 
4410  //Now add the parameter
4411  problem_pt->Dof_pt.push_back(parameter_pt);
4412  //Finally add the frequency
4413  problem_pt->Dof_pt.push_back(&Omega);
4414 
4415  // rebuild the Dof_distribution_pt
4417  Ndof*3+2,false);
4418  //Remove all previous sparse storage used during Jacobian assembly
4420  }
4421 
4422 
4423 
4424  //=======================================================================
4425  /// Destructor return the problem to its original (non-augmented) state
4426  //=======================================================================
4428  {
4429  //If we are using the block solver reset the problem's linear solver
4430  //to the original one
4431  BlockHopfLinearSolver* block_hopf_solver_pt =
4433  if(block_hopf_solver_pt)
4434  {
4435  //Reset the problem's linear solver
4436  Problem_pt->linear_solver_pt() = block_hopf_solver_pt->linear_solver_pt();
4437  //Delete the block solver
4438  delete block_hopf_solver_pt;
4439  }
4440  //Now return the problem to its original size
4441  Problem_pt->Dof_pt.resize(Ndof);
4443  Ndof,false);
4444  //Remove all previous sparse storage used during Jacobian assembly
4446  }
4447 
4448 
4449  //=============================================================
4450  ///Get the number of elemental degrees of freedom
4451  //=============================================================
4452  unsigned HopfHandler::ndof(GeneralisedElement* const &elem_pt)
4453  {
4454  unsigned raw_ndof = elem_pt->ndof();
4455  switch(Solve_which_system)
4456  {
4457  //Full augmented system
4458  case 0:
4459  return (3*raw_ndof + 2);
4460  break;
4461  //Standard non-augmented system
4462  case 1:
4463  return raw_ndof;
4464  break;
4465  //Complex system
4466  case 2:
4467  return (2*raw_ndof);
4468  break;
4469 
4470  default:
4471  throw OomphLibError("Solve_which_system can only be 0,1 or 2",
4472  OOMPH_CURRENT_FUNCTION,
4473  OOMPH_EXCEPTION_LOCATION);
4474  }
4475  }
4476 
4477  //=============================================================
4478  ///Get the global equation number of the local unknown
4479  //============================================================
4480  unsigned long HopfHandler::eqn_number(GeneralisedElement* const &elem_pt,
4481  const unsigned &ieqn_local)
4482  {
4483  //Get the raw value
4484  unsigned raw_ndof = elem_pt->ndof();
4485  unsigned long global_eqn;
4486  if(ieqn_local < raw_ndof)
4487  {
4488  global_eqn = elem_pt->eqn_number(ieqn_local);
4489  }
4490  else if(ieqn_local < 2*raw_ndof)
4491  {
4492  global_eqn = Ndof + elem_pt->eqn_number(ieqn_local - raw_ndof);
4493  }
4494  else if(ieqn_local < 3*raw_ndof)
4495  {
4496  global_eqn = 2*Ndof + elem_pt->eqn_number(ieqn_local - 2*raw_ndof);
4497  }
4498  else if(ieqn_local == 3*raw_ndof)
4499  {
4500  global_eqn = 3*Ndof;
4501  }
4502  else
4503  {
4504  global_eqn = 3*Ndof+1;
4505  }
4506  return global_eqn;
4507  }
4508 
4509  //==================================================================
4510  ///Get the residuals
4511  //=================================================================
4513  Vector<double> &residuals)
4514  {
4515  //Should only call get residuals for the full system
4516  if(Solve_which_system==0)
4517  {
4518  //Need to get raw residuals and jacobian
4519  unsigned raw_ndof = elem_pt->ndof();
4520 
4521  DenseMatrix<double> jacobian(raw_ndof), M(raw_ndof);
4522  //Get the basic residuals, jacobian and mass matrix
4523  elem_pt->get_jacobian_and_mass_matrix(residuals,jacobian,M);
4524 
4525  //Initialise the pen-ultimate residual
4526  residuals[3*raw_ndof] = -1.0/
4527  (double)(Problem_pt->mesh_pt()->nelement());
4528  residuals[3*raw_ndof+1] = 0.0;
4529 
4530  //Now multiply to fill in the residuals
4531  for(unsigned i=0;i<raw_ndof;i++)
4532  {
4533  residuals[raw_ndof+i] = 0.0;
4534  residuals[2*raw_ndof+i] = 0.0;
4535  for(unsigned j=0;j<raw_ndof;j++)
4536  {
4537  unsigned global_unknown = elem_pt->eqn_number(j);
4538  //Real part
4539  residuals[raw_ndof+i] +=
4540  jacobian(i,j)*Phi[global_unknown] + Omega*M(i,j)*Psi[global_unknown];
4541  //Imaginary part
4542  residuals[2*raw_ndof+i] +=
4543  jacobian(i,j)*Psi[global_unknown] - Omega*M(i,j)*Phi[global_unknown];
4544  }
4545  //Get the global equation number
4546  unsigned global_eqn = elem_pt->eqn_number(i);
4547 
4548  //Real part
4549  residuals[3*raw_ndof] += (Phi[global_eqn]*C[global_eqn])/
4550  Count[global_eqn];
4551  //Imaginary part
4552  residuals[3*raw_ndof+1] += (Psi[global_eqn]*C[global_eqn])/
4553  Count[global_eqn];
4554  }
4555  }
4556  else
4557  {
4558  throw OomphLibError("Solve_which_system can only be 0",
4559 OOMPH_CURRENT_FUNCTION,
4560  OOMPH_EXCEPTION_LOCATION);
4561 
4562  }
4563  }
4564 
4565 
4566  //===============================================================
4567  /// \short Calculate the elemental Jacobian matrix "d equation
4568  /// / d variable".
4569 //==================================================================
4571  Vector<double> &residuals,
4572  DenseMatrix<double> &jacobian)
4573  {
4574  //The standard case
4575  if(Solve_which_system==0)
4576  {
4577  unsigned augmented_ndof = ndof(elem_pt);
4578  unsigned raw_ndof = elem_pt->ndof();
4579 
4580  //Get the basic residuals and jacobian
4581  DenseMatrix<double> M(raw_ndof);
4582  elem_pt->get_jacobian_and_mass_matrix(residuals,jacobian,M);
4583  //Now fill in the actual residuals
4584  get_residuals(elem_pt,residuals);
4585 
4586  //Now the jacobian appears in other entries
4587  for(unsigned n=0;n<raw_ndof;++n)
4588  {
4589  for(unsigned m=0;m<raw_ndof;++m)
4590  {
4591  jacobian(raw_ndof+n,raw_ndof+m) = jacobian(n,m);
4592  jacobian(raw_ndof+n,2*raw_ndof+m) = Omega*M(n,m);
4593  jacobian(2*raw_ndof+n,2*raw_ndof+m) = jacobian(n,m);
4594  jacobian(2*raw_ndof+n,raw_ndof+m) = -Omega*M(n,m);
4595  unsigned global_eqn = elem_pt->eqn_number(m);
4596  jacobian(raw_ndof+n,3*raw_ndof+1) += M(n,m)*Psi[global_eqn];
4597  jacobian(2*raw_ndof+n,3*raw_ndof+1) -= M(n,m)*Phi[global_eqn];
4598  }
4599 
4600  unsigned local_eqn = elem_pt->eqn_number(n);
4601  jacobian(3*raw_ndof,raw_ndof+n) = C[local_eqn]/Count[local_eqn];
4602  jacobian(3*raw_ndof+1,2*raw_ndof+n) = C[local_eqn]/Count[local_eqn];
4603  }
4604 
4605  const double FD_step = 1.0e-8;
4606 
4607  Vector<double> newres_p(augmented_ndof), newres_m(augmented_ndof);
4608 
4609  //Loop over the dofs
4610  for(unsigned n=0;n<raw_ndof;n++)
4611  {
4612  //Just do the x's
4613  unsigned long global_eqn = eqn_number(elem_pt,n);
4614  double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
4615  double init = *unknown_pt;
4616  *unknown_pt += FD_step;
4617 
4618  //Get the new residuals
4619  get_residuals(elem_pt,newres_p);
4620 
4621  //Reset
4622  *unknown_pt = init;
4623 
4624  //Subtract
4625  *unknown_pt -= FD_step;
4626  get_residuals(elem_pt,newres_m);
4627 
4628  for(unsigned m=0;m<raw_ndof;m++)
4629  {
4630  jacobian(raw_ndof+m,n) =
4631  (newres_p[raw_ndof+m] - residuals[raw_ndof+m])/(FD_step);
4632  jacobian(2*raw_ndof+m,n) =
4633  (newres_p[2*raw_ndof+m] - residuals[2*raw_ndof+m])/(FD_step);
4634  }
4635  //Reset the unknown
4636  *unknown_pt = init;
4637  }
4638 
4639  {
4640  //Now do the global parameter
4641  double* unknown_pt = Problem_pt->Dof_pt[3*Ndof];
4642  double init = *unknown_pt;
4643  *unknown_pt += FD_step;
4644 
4646  //Get the new residuals
4647  get_residuals(elem_pt,newres_p);
4648 
4649  //Reset
4650  *unknown_pt = init;
4651 
4652  //Subtract
4653  *unknown_pt -= FD_step;
4654  get_residuals(elem_pt,newres_m);
4655 
4656  for(unsigned m=0;m<augmented_ndof-2;m++)
4657  {
4658  jacobian(m,3*raw_ndof) =
4659  (newres_p[m] - residuals[m])/FD_step;
4660  }
4661  //Reset the unknown
4662  *unknown_pt = init;
4664  }
4665  } //End of standard case
4666  //Normal case
4667  else if(Solve_which_system==1)
4668  {
4669  //Just get the normal jacobian and residuals
4670  elem_pt->get_jacobian(residuals,jacobian);
4671  }
4672  //Otherwise the augmented complex case
4673  else if(Solve_which_system==2)
4674  {
4675  unsigned raw_ndof = elem_pt->ndof();
4676 
4677  //Get the basic residuals and jacobian
4678  DenseMatrix<double> M(raw_ndof);
4679  elem_pt->get_jacobian_and_mass_matrix(residuals,jacobian,M);
4680 
4681  //We now need to fill in the other blocks
4682  for(unsigned n=0;n<raw_ndof;n++)
4683  {
4684  for(unsigned m=0;m<raw_ndof;m++)
4685  {
4686  jacobian(n,raw_ndof+m) = Omega*M(n,m);
4687  jacobian(raw_ndof+n,m) = -Omega*M(n,m);
4688  jacobian(raw_ndof+n,raw_ndof+m) = jacobian(n,m);
4689  }
4690  }
4691 
4692  //Now overwrite to fill in the residuals
4693  //The decision take is to solve for the mass matrix multiplied
4694  //terms in the residuals because they require no additional
4695  //information to assemble.
4696  for(unsigned n=0;n<raw_ndof;n++)
4697  {
4698  residuals[n] = 0.0;
4699  residuals[raw_ndof+n] = 0.0;
4700  for(unsigned m=0;m<raw_ndof;m++)
4701  {
4702  unsigned global_unknown = elem_pt->eqn_number(m);
4703  //Real part
4704  residuals[n] += M(n,m)*Psi[global_unknown];
4705  //Imaginary part
4706  residuals[raw_ndof+n] -= M(n,m)*Phi[global_unknown];
4707  }
4708  }
4709  } //End of complex augmented case
4710  else
4711  {
4712  throw OomphLibError("Solve_which_system can only be 0,1 or 2",
4713  OOMPH_CURRENT_FUNCTION,
4714  OOMPH_EXCEPTION_LOCATION);
4715  }
4716  }
4717 
4718 
4719 
4720  //==================================================================
4721  ///Get the derivatives of the augmented residuals with respect to
4722  ///a parameter
4723  //=================================================================
4725  GeneralisedElement* const &elem_pt,
4726  double* const &parameter_pt, Vector<double> &dres_dparam)
4727  {
4728  //Should only call get residuals for the full system
4729  if(Solve_which_system==0)
4730  {
4731  //Need to get raw residuals and jacobian
4732  unsigned raw_ndof = elem_pt->ndof();
4733 
4734  DenseMatrix<double> djac_dparam(raw_ndof), dM_dparam(raw_ndof);
4735  //Get the basic residuals, jacobian and mass matrix
4737  parameter_pt,dres_dparam,djac_dparam,dM_dparam);
4738 
4739  //Initialise the pen-ultimate residual, which does not
4740  //depend on the parameter
4741  dres_dparam[3*raw_ndof] = 0.0;
4742  dres_dparam[3*raw_ndof+1] = 0.0;
4743 
4744  //Now multiply to fill in the residuals
4745  for(unsigned i=0;i<raw_ndof;i++)
4746  {
4747  dres_dparam[raw_ndof+i] = 0.0;
4748  dres_dparam[2*raw_ndof+i] = 0.0;
4749  for(unsigned j=0;j<raw_ndof;j++)
4750  {
4751  unsigned global_unknown = elem_pt->eqn_number(j);
4752  //Real part
4753  dres_dparam[raw_ndof+i] +=
4754  djac_dparam(i,j)*Phi[global_unknown] +
4755  Omega*dM_dparam(i,j)*Psi[global_unknown];
4756  //Imaginary part
4757  dres_dparam[2*raw_ndof+i] +=
4758  djac_dparam(i,j)*Psi[global_unknown] -
4759  Omega*dM_dparam(i,j)*Phi[global_unknown];
4760  }
4761  }
4762  }
4763  else
4764  {
4765  throw OomphLibError("Solve_which_system can only be 0",
4766 OOMPH_CURRENT_FUNCTION,
4767  OOMPH_EXCEPTION_LOCATION);
4768 
4769  }
4770  }
4771 
4772 
4773 
4774 
4775  //========================================================================
4776  /// Overload the derivative of the residuals and jacobian
4777  /// with respect to a parameter so that it breaks because it should not
4778  /// be required
4779  //========================================================================
4781  GeneralisedElement* const &elem_pt,
4782  double* const &parameter_pt,
4783  Vector<double> &dres_dparam,
4784  DenseMatrix<double> &djac_dparam)
4785  {
4786  std::ostringstream error_stream;
4787  error_stream <<
4788  "This function has not been implemented because it is not required\n";
4789  error_stream << "in standard problems.\n";
4790  error_stream <<
4791  "If you find that you need it, you will have to implement it!\n\n";
4792 
4793  throw OomphLibError(error_stream.str(),
4794 OOMPH_CURRENT_FUNCTION,
4795  OOMPH_EXCEPTION_LOCATION);
4796  }
4797 
4798 
4799  //=====================================================================
4800  /// Overload the hessian vector product function so that
4801  /// it breaks because it should not be required
4802  //========================================================================
4804  GeneralisedElement* const &elem_pt,
4805  Vector<double> const &Y,
4806  DenseMatrix<double> const &C,
4807  DenseMatrix<double> &product)
4808  {
4809  std::ostringstream error_stream;
4810  error_stream <<
4811  "This function has not been implemented because it is not required\n";
4812  error_stream << "in standard problems.\n";
4813  error_stream <<
4814  "If you find that you need it, you will have to implement it!\n\n";
4815 
4816  throw OomphLibError(error_stream.str(),
4817 OOMPH_CURRENT_FUNCTION,
4818  OOMPH_EXCEPTION_LOCATION);
4819  }
4820 
4821 
4822  //==========================================================================
4823  /// Return the eigenfunction(s) associated with the bifurcation that
4824  /// has been detected in bifurcation tracking problems
4825  //==========================================================================
4827  Vector<DoubleVector> &eigenfunction)
4828  {
4829  //There is a real and imaginary part of the null vector
4830  eigenfunction.resize(2);
4832  //Rebuild the vector
4833  eigenfunction[0].build(&dist,0.0);
4834  eigenfunction[1].build(&dist,0.0);
4835  //Set the value to be the null vector
4836  for(unsigned n=0;n<Ndof;n++)
4837  {
4838  eigenfunction[0][n] = Phi[n];
4839  eigenfunction[1][n] = Psi[n];
4840  }
4841  }
4842 
4843 
4844  //====================================================================
4845  /// Set to solve the standard (underlying jacobian) system
4846  //===================================================================
4848  {
4849  if(Solve_which_system!=1)
4850  {
4851  Solve_which_system = 1;
4852  //Restrict the problem to the standard variables only
4853  Problem_pt->Dof_pt.resize(Ndof);
4855  Ndof,false);
4856  //Remove all previous sparse storage used during Jacobian assembly
4858  }
4859  }
4860 
4861  //====================================================================
4862  /// Set to solve the complex (jacobian and mass matrix) system
4863  //===================================================================
4865  {
4866  //If we were not solving the complex system resize the unknowns
4867  //accordingly
4868  if(Solve_which_system!=2)
4869  {
4871  //Resize to the first Ndofs (will work whichever system we were
4872  //solving before)
4873  Problem_pt->Dof_pt.resize(Ndof);
4874  //Add the first (real) part of the eigenfunction back into the problem
4875  for(unsigned n=0;n<Ndof;n++)
4876  {
4877  Problem_pt->Dof_pt.push_back(&Phi[n]);
4878  }
4880  Ndof*2,false);
4881  //Remove all previous sparse storage used during Jacobian assembly
4883  }
4884  }
4885 
4886 
4887  //=================================================================
4888  /// Set to Solve full system system
4889  //=================================================================
4891  {
4892  //If we are starting from another system
4893  if(Solve_which_system)
4894  {
4895  Solve_which_system = 0;
4896 
4897  //Resize to the first Ndofs (will work whichever system we were
4898  //solving before)
4899  Problem_pt->Dof_pt.resize(Ndof);
4900  //Add the additional unknowns back into the problem
4901  for(unsigned n=0;n<Ndof;n++)
4902  {
4903  Problem_pt->Dof_pt.push_back(&Phi[n]);
4904  }
4905  for(unsigned n=0;n<Ndof;n++)
4906  {
4907  Problem_pt->Dof_pt.push_back(&Psi[n]);
4908  }
4909  //Now add the parameter
4910  Problem_pt->Dof_pt.push_back(Parameter_pt);
4911  //Finally add the frequency
4912  Problem_pt->Dof_pt.push_back(&Omega);
4913 
4914  //
4916  3*Ndof+2,false);
4917  //Remove all previous sparse storage used during Jacobian assembly
4919  }
4920  }
4921 }
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.
virtual void get_residuals(Vector< double > &residuals)
Calculate the vector of residuals of the equations in the element. By default initialise the vector t...
Definition: elements.h:975
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
virtual void solve(Problem *const &problem_pt, DoubleVector &result)=0
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
bool built() const
unsigned nrow_local() const
access function for the num of local rows on this processor.
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.
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero...
double * global_dof_pt(const unsigned &i)
Return a pointer to the dof, indexed by global equation number which may be haloed or stored locally...
Definition: problem.h:1648
void solve_block_system()
Set to solve the block system.
void get_hessian_vector_products(DoubleVectorWithHaloEntries const &Y, Vector< DoubleVectorWithHaloEntries > const &C, Vector< DoubleVectorWithHaloEntries > &product)
Return the product of the global hessian (derivative of Jacobian matrix with respect to all variables...
Definition: problem.cc:7967
~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...
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the vector are redistributed to match the new distribution. In a non-MPI rebuild this...
virtual void get_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable".
Definition: elements.h:984
virtual void actions_after_change_in_bifurcation_parameter()
Actions that are to be performed after a change in the parameter that is being varied as part of the ...
Definition: problem.h:1134
virtual void get_djacobian_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Calculate the derivatives of the elemental Jacobian matrix and residuals with respect to a parameter...
Definition: elements.h:1040
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
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.
virtual void get_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
Calculate the residuals and the elemental "mass" matrix, the matrix that multiplies the time derivati...
Definition: elements.h:997
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Get the global equation number of the local unknown.
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned...
bool Enable_resolve
Boolean that indicates whether the matrix (or its factors, in the case of direct solver) should be st...
Definition: linear_solver.h:80
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.
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
LinearAlgebraDistribution * Dof_distribution_pt
The distribution of the DOFs in this problem. This object is created in the Problem constructor and s...
Definition: problem.h:457
virtual void resolve(const DoubleVector &rhs, DoubleVector &result)
Resolve the system defined by the last assembled jacobian and the rhs vector. Solution is returned in...
double *& dof_pt(const unsigned &i)
Pointer to i-th dof in the problem.
Definition: problem.h:1708
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.
double & global_value(const unsigned &i)
Direct access to global entry.
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 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 build(const OomphCommunicator *const comm_pt, const unsigned &first_row, const unsigned &nrow_local, const unsigned &nrow=0)
Sets the distribution. Takes first_row, nrow_local and nrow as arguments. If nrow is not provided or ...
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.
double dot(const DoubleVector &vec) const
compute the dot product of this vector with the vector vec.
bool is_resolve_enabled() const
Boolean flag indicating if resolves are enabled.
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.
bool distributed() const
distribution is serial or distributed
void get_derivative_wrt_global_parameter(double *const &parameter_pt, DoubleVector &result)
Get the derivative of the entire residuals vector wrt a global parameter, used in continuation proble...
Definition: problem.cc:7868
virtual void get_inner_products(Vector< std::pair< unsigned, unsigned > > const &history_index, Vector< double > &inner_product)
Return the vector of inner product of the given pairs of history values.
Definition: elements.h:1089
Problem * Problem_pt
Pointer to the problem, used in the resolve.
virtual void get_inner_product_vectors(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...
Definition: elements.h:1100
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1264
Vector< double > C
A constant vector used to ensure that the null vector is not trivial.
LinearAlgebraDistribution * Augmented_dof_distribution_pt
The augmented distribution.
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.
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".
e
Definition: cfortran.h:575
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1221
DoubleVector * Alpha_pt
Pointer to the storage for the vector alpha.
LinearAlgebraDistribution * Dof_distribution_pt
Store the original dof distribution.
virtual void get_hessian_vector_products(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...
Definition: elements.h:1077
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.
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.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
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 dof_pt_vector(Vector< double * > &dof_pt)
Return the vector of pointers to dof values.
Definition: elements.h:863
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...
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.
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
double norm() const
compute the 2 norm of this vector
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.
AssemblyHandler *& assembly_handler_pt()
Return a pointer to the assembly handler object.
Definition: problem.h:1502
~AugmentedBlockPitchForkLinearSolver()
Destructor: clean up the allocated memory.
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...
virtual void get_djacobian_and_dmass_matrix_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
Calculate the derivatives of the elemental Jacobian matrix mass matrix and residuals with respect to ...
Definition: elements.h:1055
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 build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
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.
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
double & dof(const unsigned &i)
i-th dof in the problem
Definition: problem.h:1696
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
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.
void build_halo_scheme(DoubleVectorHaloScheme *const &halo_scheme_pt)
Construct the halo scheme and storage for the halo data.
int & sign_of_jacobian()
Access function for the sign of the global jacobian matrix. This will be set by the linear solver...
Definition: problem.h:2301
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 ndof() const
Return the number of dofs.
Definition: problem.h:1574
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.
virtual void get_dresiduals_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam)
Calculate the derivatives of the residuals with respect to a parameter.
Definition: elements.h:1028
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.
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.
virtual void get_residuals(DoubleVector &residuals)
Return the fully-assembled residuals Vector for the problem: Virtual so it can be overloaded in for m...
Definition: problem.cc:3671
void setup_dof_halo_scheme()
Function that is used to setup the halo scheme.
Definition: problem.cc:284
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 get_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Calculate the residuals and jacobian and elemental "mass" matrix, the matrix that multiplies the time...
Definition: elements.h:1010
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.
virtual void enable_resolve()
Enable resolve (i.e. store matrix and/or LU decomposition, say) Virtual so it can be overloaded to pe...
Vector< Vector< unsigned > > Sparse_assemble_with_arrays_previous_allocation
the number of elements in each row of a compressed matrix in the previous matrix assembly.
Definition: problem.h:661
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.
HopfHandler(Problem *const &problem_pt, double *const &parameter_pt)
Constructor.
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
Problem * Problem_pt
Pointer to the problem, used in the resolve.
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.
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: problem.h:1424
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 synchronise()
Synchronise the halo data.
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
void solve_complex_system()
Set to solve the complex system.
Problem * Problem_pt
Pointer to the problem.
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.
Problem * Problem_pt
Pointer to the problem.
double * Parameter_pt
Storage for the pointer to the parameter.
void dof_vector(const unsigned &t, Vector< double > &dof)
Return the vector of dof values at time level t.
Definition: elements.h:837
~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.
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
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.
Vector< double * > Dof_pt
Vector of pointers to dofs.
Definition: problem.h:551
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.
DoubleVectorHaloScheme * Halo_scheme_pt
Pointer to the halo scheme for any global vectors that have the Dof_distribution. ...
Definition: problem.h:574
virtual void actions_before_newton_convergence_check()
Any actions that are to be performed before the residual is checked in the Newton method...
Definition: problem.h:1031
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.
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 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.
bool distributed() const
Definition: problem.h:798
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.
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...
bool is_halo() const
Is this element a halo?
Definition: elements.h:1141
Vector< double > Psi
The imaginary part of the null vector.
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.