problem.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1302 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-10-23 12:52:19 +0100 (Mon, 23 Oct 2017) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //A generic problem class to keep MH happy
31 
32 //Include guards to prevent multiple inclusion of this header
33 #ifndef OOMPH_PROBLEM_CLASS_HEADER
34 #define OOMPH_PROBLEM_CLASS_HEADER
35 
36 
37 // Config header generated by autoconfig
38 #ifdef HAVE_CONFIG_H
39  #include <oomph-lib-config.h>
40 #endif
41 
42 #ifdef OOMPH_HAS_MPI
43 #include "mpi.h"
44 #endif
45 
46 //OOMPH-LIB headers
47 #include "Vector.h"
48 #include "matrices.h"
50 #include "explicit_timesteppers.h"
52 #include <complex>
53 #include <map>
54 
55 
56 namespace oomph
57 {
58  //Forward definition for Data class
59  class Data;
60 
61  //Forward definition for Time class
62  class Time;
63 
64  //Forward definition for TimeStepper class
65  class TimeStepper;
66 
67  //Forward definition for Mesh class
68  class Mesh;
69 
70  //Forward definition for RefineableElement class
71  class RefineableElement;
72 
73  //Forward definition for PRefineableElement class
74  class PRefineableElement;
75 
76  //Forward definition for FiniteElement class
77  class FiniteElement;
78 
79  //Forward definition for the GeneralisedElement class
80  class GeneralisedElement;
81 
82  //Forward definition for the Linearsolver class
83  class LinearSolver;
84 
85  //Forward definition for the Eigensolver class
86  class EigenSolver;
87 
88  //Forward definition for the assembly handler
89  class AssemblyHandler;
90 
91  //Forward definition for sum of matrices class
92  class SumOfMatrices;
93 
94 /////////////////////////////////////////////////////////////////////
95 /////////////////////////////////////////////////////////////////////
96 /////////////////////////////////////////////////////////////////////
97 
98 
99 //=======================================================================
100 /// \short The Problem class
101 ///
102 /// The main components of a Problem are:
103 /// - a pointer to the (global) Mesh (which provides ordered
104 /// access to the nodes, elements, etc of all submeshes in the problem --
105 /// if there are any)
106 /// - pointers to the submeshes (if there are any)
107 /// - pointers to any global Data
108 /// - a pointer to the global Time
109 /// - pointers to the Timestepper used in the problem
110 /// - a pointer to the linear solver that will be used in Newton's method
111 /// - pointers to all DOFs in the problem.
112 ///
113 /// Obviously, at this level in the code hierarchy, many things become
114 /// very problem-dependent but when setting up a specific problem
115 /// (as a class that inherits from Problem), the problem constructor
116 /// will/should typically have a structure similar to this:
117 /// \code
118 /// // Set up a timestepper
119 /// TimeStepper* time_stepper_pt=new Steady;
120 ///
121 /// // Add timestepper to problem
122 /// add_time_stepper_pt(time_stepper_pt);
123 ///
124 /// // Build and assign mesh
125 /// Problem::mesh_pt() = new SomeMesh<ELEMENT_TYPE>(Nelement,time_stepper_pt);
126 ///
127 /// // Set the boundary conditions for this problem: All nodes are
128 /// // free by default -- just pin the ones that have Dirichlet conditions
129 /// // here (boundary 0 in this example)
130 /// unsigned i_bound = 0
131 /// unsigned num_nod= mesh_pt()->nboundary_node(ibound);
132 /// for (unsigned inod=0;inod<num_nod;inod++)
133 /// {
134 /// mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
135 /// }
136 ///
137 /// // Complete the build of all elements so they are fully functional
138 ///
139 /// // Setup equation numbering scheme
140 /// oomph_info <<"Number of equations: " << assign_eqn_numbers() << std::endl;
141 /// \endcode
142 /// For time-dependent problems, we can then use
143 /// \code assign_initial_values_impulsive (); \endcode
144 /// to generate an initial guess (for an impulsive start) and
145 /// the problem can then be solved, e.g. using the steady
146 /// or unsteady Newton solvers
147 /// \code newton_solve() \endcode
148 /// or
149 /// \code unsteady_newton_solve(...) \endcode
150 ///
151 //=======================================================================
153 {
154  //The handler classes need to be friend
155  friend class FoldHandler;
156  friend class PitchForkHandler;
157  friend class HopfHandler;
158  template<unsigned NNODE_1D>
160  friend class BlockFoldLinearSolver;
164  friend class BlockHopfLinearSolver;
165 
166 
167  private:
168 
169  /// The mesh pointer
171 
172  /// Vector of pointers to submeshes
174 
175  /// Pointer to the linear solver for the problem
177 
178  /// Pointer to the linear solver used for explicit time steps (this is
179  /// likely to be different to the linear solver for newton solves because
180  /// explicit time steps only involve inverting a mass matrix. This can be
181  /// done very efficiently by, e.g. CG with a diagonal predconditioner).
183 
184  /// Pointer to the eigen solver for the problem
186 
187  //Pointer to handler
189 
190  /// Pointer to the default linear solver
192 
193  /// Pointer to the default eigensolver
195 
196  /// Pointer to the default assembly handler
198 
199  /// Pointer to global time for the problem
201 
202  /// \short The Vector of time steppers (there could be many
203  /// different ones in multiphysics problems)
205 
206  /// \short Pointer to a single explicit timestepper
208 
209  /// Pointer to vector for backup of dofs
211 
212  /// \short Has default set_initial_condition function been called?
213  /// Default: false
215 
216  /// \short Use the globally convergent newton method
218 
219  /// \short Boolean to indicate that empty
220  /// actions_before_read_unstructured_meshes() function has been called.
222 
223  /// \short Boolean to indicate that empty
224  /// actions_after_read_unstructured_meshes() function has been called.
226 
227  /// \short Boolean to indicate whether local dof pointers should be
228  /// stored in the elements
230 
231  /// \short Use values from the time stepper predictor as an initial guess
233 
234  protected:
235 
236  ///\short Vector of pointers to copies of the problem used in adaptive
237  ///bifurcation tracking problems (ALH: TEMPORARY HACK, WILL BE FIXED)
239 
240  /// \short Map used to determine whether the derivatives with respect to
241  /// a parameter should be finite differenced. The default is that
242  /// finite differences should be used
243  std::map<double*,bool> Calculate_dparameter_analytic;
244 
245  /// \short Map used to determine whether the hessian products should be
246  /// computed using finite differences. The default is that finite differences
247  /// will be used
249 
250  public:
251 
252  /// \short Hook for debugging. Can be overloaded in driver code; argument
253  /// allows identification of where we're coming from
254  virtual void debug_hook_fct(const unsigned& i)
255  {
256  oomph_info << "Called empty hook fct with i=" << i << std::endl;
257  }
258 
259  /// \short Function to turn on analytic calculation of the parameter
260  /// derivatives in continuation and bifurcation detection problems
261  inline void set_analytic_dparameter(double* const &parameter_pt)
262  {Calculate_dparameter_analytic[parameter_pt] = true;}
263 
264  /// \short Function to turn off analytic calculation of the parameter
265  /// derivatives in continuation and bifurcation detection problems
266  inline void unset_analytic_dparameter(double* const &parameter_pt)
267  {
268  //Find the iterator to the parameter
269  std::map<double*,bool>::iterator it =
270  Calculate_dparameter_analytic.find(parameter_pt);
271  //If the parameter has been found, erase it
272  if(it!=Calculate_dparameter_analytic.end())
273  {
275  }
276  }
277 
278  /// \short Function to determine whether the parameter derivatives
279  /// are calculated analytically
280  inline bool is_dparameter_calculated_analytically(double* const &parameter_pt)
281  {
282  //If the entry is found then the iterator returned from the find
283  //command will NOT be equal to the end and the expression will
284  //return true, otherwise it will return false
285  return (Calculate_dparameter_analytic.find(parameter_pt)!=
287  }
288 
289  /// \short Function to turn on analytic calculation of the parameter
290  /// derivatives in continuation and bifurcation detection problems
293 
294  /// \short Function to turn off analytic calculation of the parameter
295  /// derivatives in continuation and bifurcation detection problems
298 
299  /// \short Function to determine whether the hessian products
300  /// are calculated analytically
303 
304  /// \short Set all pinned values to zero.
305  /// Used to set boundary conditions to be homogeneous in the copy
306  /// of the problem used in adaptive bifurcation tracking
307  /// (ALH: TEMPORARY HACK, WILL BE FIXED)
309 
310 
311  /// \short Flag to allow suppression of warning messages re reading in
312  /// unstructured meshes during restart.
314 
315 
316  private:
317 
318 
319  /// \short Private helper function that actually performs the unsteady
320  /// "doubly" adaptive Newton solve. See actual (non-helper) functions
321  /// for description of parameters.
323  const double &dt,
324  const double &epsilon,
325  const unsigned &max_adapt,
326  const unsigned& suppress_resolve_after_spatial_adapt,
327  const bool &first,
328  const bool& shift=true);
329 
330 
331  /// \short Helper function to do compund refinement of (all) refineable
332  /// (sub)mesh(es) uniformly as many times as specified in vector and
333  /// rebuild problem; doc refinement process. Set boolean argument
334  /// to true if you want to prune immediately after refining the meshes
335  /// individually.
336  void refine_uniformly_aux(const Vector<unsigned>& nrefine_for_mesh,
337  DocInfo& doc_info,
338  const bool& prune);
339 
340  /// \short Helper function to do compund p-refinement of (all) p-refineable
341  /// (sub)mesh(es) uniformly as many times as specified in vector and
342  /// rebuild problem; doc refinement process. Set boolean argument
343  /// to true if you want to prune immediately after refining the meshes
344  /// individually.
345  void p_refine_uniformly_aux(const Vector<unsigned>& nrefine_for_mesh,
346  DocInfo& doc_info,
347  const bool& prune);
348 
349  /// \short Helper function to re-setup the Base_mesh enumeration
350  /// (used during load balancing) after pruning
352 
353  /// \short Private helper function that is used to assemble the Jacobian
354  /// matrix in the case when the storage is row or column compressed.
355  /// The boolean Flag indicates
356  /// if we want compressed row format (true) or compressed column.
357  /// This version uses vectors of pairs.
359  Vector<int* > &column_or_row_index,
360  Vector<int* > &row_or_column_start,
361  Vector<double* > &value,
362  Vector<unsigned> &nnz,
363  Vector<double* > &residual,
364  bool compressed_row_flag);
365 
366  /// \short Private helper function that is used to assemble the Jacobian
367  /// matrix in the case when the storage is row or column compressed.
368  /// The boolean Flag indicates
369  /// if we want compressed row format (true) or compressed column.
370  /// This version uses two vectors.
372  Vector<int* > &column_or_row_index,
373  Vector<int* > &row_or_column_start,
374  Vector<double* > &value,
375  Vector<unsigned> &nnz,
376  Vector<double* > &residual,
377  bool compressed_row_flag);
378 
379  /// \short Private helper function that is used to assemble the Jacobian
380  /// matrix in the case when the storage is row or column compressed.
381  /// The boolean Flag indicates
382  /// if we want compressed row format (true) or compressed column.
383  /// This version uses maps
385  Vector<int* > &column_or_row_index,
386  Vector<int* > &row_or_column_start,
387  Vector<double* > &value,
388  Vector<unsigned> &nnz,
389  Vector<double* > &residual,
390  bool compressed_row_flag);
391 
392  /// \short Private helper function that is used to assemble the Jacobian
393  /// matrix in the case when the storage is row or column compressed.
394  /// The boolean Flag indicates
395  /// if we want compressed row format (true) or compressed column.
396  /// This version uses lists
398  Vector<int* > &column_or_row_index,
399  Vector<int* > &row_or_column_start,
400  Vector<double* > &value,
401  Vector<unsigned> &nnz,
402  Vector<double* > &residual,
403  bool compressed_row_flag);
404 
405  /// \short Private helper function that is used to assemble the Jacobian
406  /// matrix in the case when the storage is row or column compressed.
407  /// The boolean Flag indicates
408  /// if we want compressed row format (true) or compressed column.
409  /// This version uses lists
411  Vector<int* > &column_or_row_index,
412  Vector<int* > &row_or_column_start,
413  Vector<double* > &value,
414  Vector<unsigned> &nnz,
415  Vector<double* > &residual,
416  bool compressed_row_flag);
417 
418  /// \short Vector of global data: "Nobody" (i.e. none of the elements etc.)
419  /// is "in charge" of this Data so it would be overlooked when it
420  /// comes to equation-numbering, timestepping etc. Including
421  /// (pointers) to such Data in here, puts the Problem itself
422  /// "in charge" of these tasks.
424 
425  /// \short A function that performs the guts of the continuation
426  /// derivative calculation in arc length continuation problems.
428 
429  /// \short A function that performs the guts of the continuation
430  /// derivative calculation in arc-length continuation problems using
431  /// finite differences
433  double* const &parameter_pt);
434 
435  /// \short A function that is used to adapt a bifurcation-tracking
436  /// problem, which requires separate interpolation of the
437  /// associated eigenfunction. The error measure is chosen to be
438  /// a suitable combination of the errors in the base flow and the
439  /// eigenfunction
440  void bifurcation_adapt_helper(unsigned &n_refined, unsigned &n_unrefined,
441  const unsigned &bifurcation_type,
442  const bool &actually_adapt=true);
443 
444  /// \short A function that is used to document the errors used
445  /// in the adaptive solution of bifurcation problems.
446  void bifurcation_adapt_doc_errors(const unsigned &bifurcation_type);
447 
448 //ALH_TEMP_DEVELOPMENT
449  protected:
450  /// \short The distribution of the DOFs in this problem.
451  /// This object is created in the Problem constructor and setup
452  /// when assign_eqn_numbers(...) is called.
453  /// If this problem is distributed then this distribution will match
454  /// the distribution of the equation numbers.
455  /// If this problem is not distributed then this distribution will
456  /// be uniform over all processors.
458 
459  private:
460 #ifdef OOMPH_HAS_MPI
461 
462  ///\short Helper method that returns the (unique) global equations to which
463  ///the elements in the range el_lo to el_hi contribute on this
464  ///processor using the given assembly_handler
466  const unsigned &el_lo, const unsigned &el_hi,
467  Vector<unsigned> &my_eqns);
468 
469  /// \short Helper method to assemble CRDoubleMatrices from distributed
470  /// on multiple processors.
472  (const LinearAlgebraDistribution* const &dist_pt,
473  Vector<int* > &column_or_row_index,
474  Vector<int* > &row_or_column_start,
475  Vector<double* > &value,
476  Vector<unsigned > &nnz,
477  Vector<double* > &residuals);
478 
479  /// \short A private helper function to
480  /// copy the haloed equation numbers into the halo equation numbers,
481  /// either for the problem's one and only mesh or for all of its
482  /// submeshes. Bools control if we deal with data associated with
483  /// external halo/ed elements/nodes or the "normal" halo/ed ones.
484  void copy_haloed_eqn_numbers_helper(const bool& do_halos,
485  const bool& do_external_halos);
486 
487  /// \short Private helper function to remove repeated data
488  /// in external haloed elements in specified mesh. Bool is true if some data
489  /// was removed -- this usually requires re-running through certain
490  /// parts of the equation numbering procedure.
491  void remove_duplicate_data(Mesh* const &mesh_pt,
492  bool& actually_removed_some_data);
493 
494 
495  /// \short Consolidate external halo node storage by removing nulled out
496  /// pointers in external halo and haloed schemes for all meshes.
498 
499  /// \short Helper function to re-assign the first and last elements to be
500  /// assembled by each processor during parallel assembly for
501  /// non-distributed problem.
503 
504  /// \short Boolean to switch on assessment of load imbalance in parallel
505  /// assembly of distributed problem
507 
508  /// \short Flag to use "default partition" during load balance.
509  /// Should only be set to true when run in validation mode.
511 
512  /// \short First element to be assembled by given processor for
513  /// non-distributed problem (only kept up to date when default assignment
514  /// is used)
516 
517  /// \short Last element (plus one) to be assembled by given processor
518  /// for non-distributed problem (only kept up to date when default
519  /// assignment is used)
521 
522  /// \short Boolean indicating that the division of elements over processors
523  /// during the assembly process must be re-load-balanced.
524  /// (only used for non-distributed problems)
526 
527  /// \short Map which stores the correspondence between a root element and
528  /// its element number (plus one) within the global mesh at the point
529  /// when it is distributed. NB a root element in this instance is one
530  /// of the elements in the uniformly-refined mesh at the point when
531  /// Problem::distribute() is called, since these elements become roots
532  /// on each of the processors involved in the distribution. Null when
533  /// element doesn't exist following the adjustment of this when pruning.
534  std::map<GeneralisedElement*,unsigned> Base_mesh_element_number_plus_one;
535 
536 
537  /// \short Vector to store the correspondence between a root element and
538  /// its element number within the global mesh at the point when it is
539  /// distributed. NB a root element in this instance is one of the elements
540  /// in the uniformly-refined mesh at the point when Problem::distribute()
541  /// is called, since these elements become roots on each of the processors
542  /// involved in the distribution. Null when element doesn't exist
543  /// following the adjustment of this when pruning.
545 
546 #endif
547 
548  protected:
549 
550  /// Vector of pointers to dofs
552 
553  ///\short Counter that records how many elements contribute to each dof.
554  /// Used to determine the (discrete) arc-length automatically.
555  /// It really should be an integer, but is a double so that the
556  /// distribution information can be used for distributed problems
558 
559  /// \short Function that populates the Element_counter_per_dof vector
560  /// with the number of elements that contribute to each dof. For example,
561  /// with linear elements in 1D each dof contains contributions from two
562  /// elements apart from those on the boundary. Returns the total number
563  /// of elements in the problem
564  unsigned setup_element_count_per_dof();
565 
566 
567 #ifdef OOMPH_HAS_MPI
568 
569  /// Storage for assembly times (used for load balancing)
571 
572  /// \short Pointer to the halo scheme for any global vectors
573  /// that have the Dof_distribution
575 
576  /// \short Storage for the halo degrees of freedom (only required)
577  /// when accessing via the global equation number in distributed problems
579 
580  /// \short Function that is used to setup the halo scheme
581  void setup_dof_halo_scheme();
582 
583 #endif
584 
585  //--------------------- Newton solver parameters
586 
587  /// \short Relaxation fator for Newton method (only a fractional Newton
588  /// correction is applied if this is less than 1; default 1).
590 
591  /// \short The Tolerance below which the Newton Method is deemed to have
592  /// converged
594 
595  /// Maximum number of Newton iterations
597 
598  /// \short Actual number of Newton iterations taken during the most recent
599  /// iteration
601 
602  /// Maximum residuals at start and after each newton iteration
604 
605  /// \short Maximum desired residual:
606  /// if the maximum residual exceeds this value, the program will exit
608 
609  /// \short Bool to specify what to do if a Newton solve fails within a
610  /// time adaptive solve. Default (false) is to half the step and try
611  /// again. If true then crash instead.
613 
614  /// Is re-use of Jacobian in Newton iteration enabled? Default: false
616 
617  /// \short Has a Jacobian been computed (and can therefore be re-used
618  /// if required)? Default: false
620 
621  /// \short Boolean flag indicating if we're dealing with a linear or nonlinear
622  /// Problem -- if set to false the Newton solver will not check
623  /// the residual before or after the linear solve. Set to true by default;
624  /// can be over-written in specific Problem class.
626 
627  /// \short Protected boolean flag to halt program execution
628  /// during sparse assemble process to assess peak memory
629  /// usage. Initialised to false (obviously!)
631 
632  /// \short Protected boolean flag to provide comprehensive timimings
633  /// during problem distribution. Initialised to false.
635 
636  /// \short Flag to determine which sparse assembly method to use
637  /// By default we use assembly by vectors of pairs.
639 
640  /// Enumerated flags to determine which sparse assembly method is used
646 
647  /// \short the number of elements to initially allocate for a matrix row
648  /// within the sparse_assembly_with_two_arrays(...) method (if a matrix
649  /// of this size has not been assembled already). If a matrix of this size
650  /// has been assembled then the number of elements in each row in that matrix
651  /// is used as the initial allocation
653 
654  /// \short the number of elements to add to a matrix row when the initial
655  /// allocation is exceeded within the sparse_assemble_with_two_arrays(...)
656  /// method.
658 
659  /// \short the number of elements in each row of a compressed matrix
660  /// in the previous matrix assembly.
662 
663  /// \short A tolerance used to determine whether the entry in a sparse
664  /// matrix is zero. If it is then storage need not be allocated.
666 
667  /// \short Protected helper function that is used to assemble the Jacobian
668  /// matrix in the case when the storage is row or column compressed.
669  /// The boolean Flag indicates
670  /// if we want compressed row format (true) or compressed column.
672  Vector<int* > &column_or_row_index,
673  Vector<int* > &row_or_column_start,
674  Vector<double* > &value,
675  Vector<unsigned> &nnz,
676  Vector<double* > &residual,
677  bool compressed_row_flag);
678 
680 
681  //---------------------Explicit time-stepping parameters
682 
683  ///Is re-use of the mass matrix in explicit timestepping enabled Default:false
685 
686  ///\short Has the mass matrix been computed (and can therefore be reused)
687  ///Default: false
689 
690  //---------------------Discontinuous control flags
691  ///\short Is the problem a discontinuous one, i.e. can the
692  ///elemental contributions be treated independently. Default: false
694 
695 
696  //--------------------- Adaptive time-stepping parameters
697 
698  /// Minimum desired dt: if dt falls below this value, exit
699  double Minimum_dt;
700 
701  /// Maximum desired dt
702  double Maximum_dt;
703 
704  /// \short Maximum possible increase of dt between time-steps in adaptive
705  /// schemes
707 
708  /// \short Minimum allowed decrease of dt between time-steps in adaptive
709  /// schemes. Lower scaling values will reject the time-step (and retry
710  /// with a smaller dt).
712 
713  /// \short If Minimum_dt_but_still_proceed positive, then dt will not be
714  /// reduced below this value during adaptive timestepping and the
715  /// computation will continue with this value, accepting the larger
716  /// errors that this will incur). Note: This option is disabled by default
717  /// as this value is initialised to -1.0.
719 
720 
721  //--------------------- Arc-length continuation paramaters
722 
723  /// Boolean to control whether arc-length should be scaled
725 
726  /// Proportion of the arc-length to taken by the parameter
728 
729  /// \short Value of the scaling parameter required so that the parameter
730  /// occupies the desired proportion of the arc length
732 
733  /// Storage for the sign of the global Jacobian
735 
736  /// \short The direction of the change in parameter that will ensure that
737  /// a branch is followed in one direction only
739 
740  /// Storage for the derivative of the global parameter wrt arc-length
742 
743  /// Storage for the present value of the global parameter
745 
746  /// Boolean to control original or new storage of dof stuff
748 
749  /// Storage for the single static continuation timestorage object
751 
752  /// Storage for the offset for the continuation derivatives from the
753  /// original dofs (default is 1, but this will be differnet when continuing
754  /// bifurcations and periodic orbits)
756 
757  /// Storage for the offset for the current values of dofs from the
758  /// original dofs (default is 2, but this will be differnet when continuing
759  /// bifurcations and periodic orbits)
761 
762  /// Storage for the derivative of the problem variables wrt arc-length
764 
765  /// Storage for the present values of the variables
767 
768  /// Storage for the current step value
769  double Ds_current;
770 
771  /// \short The desired number of Newton Steps to reach convergence at each
772  /// step along the arc
774 
775  /// Minimum desired value of arc-length
776  double Minimum_ds;
777 
778  /// Boolean to control bifurcation detection via determinant of Jacobian
780 
781  /// Boolean to control wheter bisection is used to located bifurcation
783 
784  /// Boolean to indicate whether a sign change has occured in the Jacobian
786 
787  /// Boolean to indicate whether an arc-length step has been taken
789 
790  /// Boolean to specify which scheme to use to calculate the continuation
791  /// derivatievs
793 
794 public:
795 
796  /// If we have MPI return the "problem has been distributed" flag,
797  /// otherwise it can't be distributed so return false.
798  bool distributed() const
799  {
800 #ifdef OOMPH_HAS_MPI
802 #else
803  return false;
804 #endif
805  }
806 
807 #ifdef OOMPH_HAS_MPI
808 
809 
810  /// \short enum for distribution of distributed jacobians.
811  /// 1 - Automatic - the Problem distribution is employed, unless any
812  /// processor has number of rows equal to 110% of N/P, in which case
813  /// uniform distribution is employed.
814  /// 2 - Problem - the jacobian on processor p only contains rows that
815  /// correspond to equations that are on this processor. (minimises
816  /// communication)
817  /// 3 - Uniform - each processor holds as close to N/P matrix rows as
818  /// possible. (very well load balanced)
823 
824  /// \short accesss function to the distributed matrix distribution method
825  /// 1 - Automatic - the Problem distribution is employed, unless any
826  /// processor has number of rows equal to 110% of N/P, in which case
827  /// uniform distribution is employed.
828  /// 2 - Problem - the jacobian on processor p only contains rows that
829  /// correspond to equations that are on this processor. (minimises
830  /// communication)
831  /// 3 - Uniform - each processor holds as close to N/P matrix rows as
832  /// possible. (very well load balanced)
835  {
837  }
838 
839  /// \short Enable doc of load imbalance in parallel
840  /// assembly of distributed problem
843 
844  /// \short Disable doc of load imbalance in parallel
845  /// assembly of distributed problem
848 
849  /// \short Return vector of most-recent elemental assembly times
850  /// (used for load balancing). Zero sized if no Jacobian has been
851  /// computed since last re-assignment of equation numbers
853  {return Elemental_assembly_time;}
854 
855  /// \short Clear storage of most-recent elemental assembly times
856  /// (used for load balancing). Next load balancing operation
857  /// will be based on the number of elements associated with a tree root.
859  {
861  Elemental_assembly_time.clear();
862  }
863 
864  private:
865 
866  /// \short Load balance helper routine: Get data to be sent to other
867  /// processors during load balancing and other information about
868  /// re-distribution.
869  /// - target_domain_for_local_non_halo_element: Input, generated by METIS.
870  /// target_domain_for_local_non_halo_element[e] contains the number
871  /// of the domain [0,1,...,nproc-1] to which non-halo element e on THE
872  /// CURRENT PROCESSOR ONLY has been assigned. The order of the non-halo
873  /// elements is the same as in the Problem's mesh, with the halo
874  /// elements being skipped.
875  /// - send_n: Output, number of data to be sent to each processor
876  /// - send_data: Output, storage for all values to be sent to all processors
877  /// - send_displacement: Output, start location within send_data for data to
878  /// be sent to each processor
879  /// - max_refinement_level_overall: Output, max. refinement level of any
880  /// element
882  (const Vector<unsigned>& element_domain_on_this_proc,
883  Vector<int>& send_n,
884  Vector<double>& send_data,
885  Vector<int>& send_displacement,
886  Vector<unsigned>& old_domain_for_base_element,
887  Vector<unsigned>& new_domain_for_base_element,
888  unsigned& max_refinement_level_overall);
889 
890  /// \short Get flat-packed refinement pattern for each root element in current
891  /// mesh (labeled by unique number of root element in unrefined base mesh).
892  /// The vector stored for each root element contains the following
893  /// information:
894  /// - First entry: Number of tree nodes (not just leaves!) in refinement
895  /// tree emanating from this root [Zero if root element is not refineable]
896  /// - Loop over all refinement levels
897  /// - Loop over all tree nodes (not just leaves!)
898  /// - If associated element exists when the mesh has been refined to
899  /// this level (either because it has been refined to this level or
900  /// because it's less refined): 1
901  /// - If the element is to be refined: 1; else: 0
902  /// - else (element doesn't exist when mesh is refined to this level
903  /// (because it's more refined): 0
904  /// .
905  /// .
906  /// .
908  const Vector<unsigned>& old_domain_for_base_element,
909  const Vector<unsigned>& new_domain_for_base_element,
910  const unsigned& max_refinement_level_overall,
911  std::map<unsigned, Vector<unsigned> >& flat_packed_refinement_info_for_root);
912 
913 
914 
915  /// \short Load balance helper routine: Send data to other
916  /// processors during load balancing.
917  /// - send_n: Input, number of data to be sent to each processor
918  /// - send_data: Input, storage for all values to be sent to all processors
919  /// - send_displacement: Input, start location within send_data for data to
920  /// be sent to each processor
922  (Vector<int>& send_n, Vector<double>& send_data,
923  Vector<int>& send_displacement);
924 
925  /// Send refinement information between processors
927  Vector<unsigned>& old_domain_for_base_element,
928  Vector<unsigned>& new_domain_for_base_element,
929  const unsigned& max_refinement_level_overall,
930  std::map<unsigned, Vector<unsigned> >& refinement_info_for_root_local,
931  Vector<Vector<Vector<unsigned> > >& refinement_info_for_root_elements);
932 
933  /// \short The distributed matrix distribution method
934  /// 1 - Automatic - the Problem distribution is employed, unless any
935  /// processor has number of rows equal to 110% of N/P, in which case
936  /// uniform distribution is employed.
937  /// 2 - Problem - the jacobian on processor p only contains rows that
938  /// correspond to equations that are on this processor. (minimises
939  /// communication)
940  /// 3 - Uniform - each processor holds as close to N/P matrix rows as
941  /// possible. (very well load balanced)
943 
944  /// \short The amount of data allocated during the previous parallel sparse
945  /// assemble. Used to optimise the next call to parallel_sparse_assemble()
947 
948  protected:
949 
950  /// Has the problem been distributed amongst multiple processors?
952 
953  /// Boolean to bypass check of increase in dofs during pruning
955 
956  public:
957 
958  /// Enable problem distributed
960 
961  /// Disable problem distributed
963 
964  /// \short Set default first and last elements for parallel assembly
965  /// of non-distributed problem.
967 
968  /// \short Manually set first and last elements for parallel assembly
969  /// of non-distributed problem.
971  Vector<unsigned>& first_el_for_assembly,
972  Vector<unsigned>& last_el_for_assembly)
973  {
974  First_el_for_assembly=first_el_for_assembly;
975  Last_el_plus_one_for_assembly=last_el_for_assembly;
976  unsigned n=last_el_for_assembly.size();
977  for (unsigned i=0;i<n;i++)
978  {
980  }
981  }
982 
983 
984  /// \short Get first and last elements for parallel assembly
985  /// of non-distributed problem.
987  Vector<unsigned>& first_el_for_assembly,
988  Vector<unsigned>& last_el_for_assembly) const
989  {
990  first_el_for_assembly=First_el_for_assembly;
991  last_el_for_assembly=Last_el_plus_one_for_assembly;
992  unsigned n=last_el_for_assembly.size();
993  for (unsigned i=0;i<n;i++)
994  {
995  last_el_for_assembly[i]-=1;
996  }
997  }
998 
999 #endif
1000 
1001  /// \short Actions that are to be performed before a mesh adaptation.
1002  /// These might include removing any additional elements, such as traction
1003  /// boundary elements before the adaptation.
1004  virtual void actions_before_adapt(){}
1005 
1006  /// Actions that are to be performed after a mesh adaptation.
1007  virtual void actions_after_adapt(){}
1008 
1009  protected:
1010 
1011  /// \short Any actions that are to be performed before a complete
1012  /// Newton solve (e.g. adjust boundary conditions). CAREFUL: This
1013  /// step should (and if the FD-based LinearSolver FD_LU is used, must)
1014  /// only update values that are pinned!
1015  virtual void actions_before_newton_solve() { }
1016 
1017  /// \short Any actions that are to be performed after a complete Newton solve,
1018  /// e.g. post processing. CAREFUL: This
1019  /// step should (and if the FD-based LinearSolver FD_LU is used, must)
1020  /// only update values that are pinned!
1021  virtual void actions_after_newton_solve() { }
1022 
1023  /// \short Any actions that are to be performed before
1024  /// the residual is checked in the Newton method, e.g. update
1025  /// any boundary conditions that depend upon
1026  /// variables of the problem; update any `slave' variables; or
1027  /// perform an update of the nodal positions in SpineMeshes etc.
1028  /// CAREFUL: This
1029  /// step should (and if the FD-based LinearSolver FD_LU is used, must)
1030  /// only update values that are pinned!
1032 
1033  /// \short Any actions that are to be performed before each individual
1034  /// Newton step. Most likely to be used for diagnostic purposes
1035  /// to doc the solution during a non-converging iteration, say.
1036  virtual void actions_before_newton_step() {}
1037 
1038  /// \short Any actions that are to be performed after each individual
1039  /// Newton step. Most likely to be used for diagnostic purposes
1040  /// to doc the solution during a non-converging iteration, say.
1041  virtual void actions_after_newton_step(){}
1042 
1043  /// \short Actions that should be performed before each implicit time step.
1044  /// This is needed when one wants
1045  /// to solve a steady problem before timestepping and needs to distinguish
1046  /// between the two cases
1048 
1049  /// \short Actions that should be performed after each implicit time step.
1050  /// This is needed when one wants
1051  /// to solve a steady problem before timestepping and needs to distinguish
1052  /// between the two cases
1054 
1055  /// \short Actions that should be performed after each implicit time step.
1056  /// This is needed if your actions_after_implicit_timestep() modify the
1057  /// solution in a way that affects the error estimate.
1059 
1060  /// \short Actions that should be performed before each explicit time step.
1062 
1063  /// \short Actions that should be performed after each explicit time step.
1065 
1066  /// \short Actions that are to be performed before reading in
1067  /// restart data for problems involving unstructured bulk meshes.
1068  /// If the problem contains such meshes we need
1069  /// to strip out any face elements that are attached to them
1070  /// because restart of unstructured meshes re-creates their elements
1071  /// and nodes from scratch, leading to dangling pointers from the
1072  /// face elements to the old elements and nodes. This function is
1073  /// virtual and (practically) empty
1074  /// but toggles a flag to indicate that it has been called. This is used to
1075  /// issue a warning, prompting the user to consider overloading it
1076  /// if the problem is found to contain unstructured bulk meshes during
1077  /// restarts.
1079  {
1081  }
1082 
1083  /// \short Actions that are to be performed before reading in
1084  /// restart data for problems involving unstructured bulk meshes.
1085  /// Typically used to re-attach FaceElements, say, that were stripped
1086  /// out in actions_before_read_unstructured_meshes(). This function is
1087  /// virtual and (practically) empty
1088  /// but toggles a flag to indicate that it has been called. This is used to
1089  /// issue a warning, prompting the user to consider overloading it
1090  /// if the problem is found to contain unstructured bulk meshes during
1091  /// restarts.
1093  {
1095  }
1096 
1097 #ifdef OOMPH_HAS_MPI
1098  /// Actions to be performed before a (mesh) distribution
1099  virtual void actions_before_distribute() {}
1100 
1101  /// Actions to be performed after a (mesh) distribution
1102  virtual void actions_after_distribute() {}
1103 
1104 #endif
1105 
1106  /// \short Actions that are to be performed when the global parameter
1107  /// addressed by parameter_pt
1108  /// has been changed in the function get_derivative_wrt_global_parameter()
1109  /// The default is to call actions_before_newton_solve(),
1110  /// actions_before_newton_convergence_check() and
1111  /// actions_after_newton_solve().
1112  /// This could be amazingly inefficient in certain problems and should
1113  /// be overloaded in such cases. An example would be when a remesh is
1114  /// required in general, but the global parameter does not affect the mesh
1115  /// directly.
1117  double* const &parameter_pt)
1118  {
1119  //Default action should cover all possibilities
1123  }
1124 
1125  /// \short Actions that are to be performed after a change in the parameter
1126  /// that is being varied as part of the solution of a bifurcation detection
1127  /// problem. The default is to call actions_before_newton_solve(),
1128  /// actions_before_newton_convergence_check() and
1129  /// actions_after_newton_solve().
1130  /// This could be amazingly inefficient in certain problems and should
1131  /// be overloaded in such cases. An example would be when a remesh is
1132  /// required in general, but the global parameter does not affect the mesh
1133  /// directly.
1135  {
1136  //Default action should cover all possibilities
1140  }
1141 
1142  /// \short Empty virtual function; provides hook to perform actions after the
1143  /// increase in the arclength parameter (during continuation)
1144  virtual void actions_after_parameter_increase(double* const &parameter_pt)
1145  {}
1146 
1147  /// \short Access function to the derivative of the i-th (local)
1148  /// dof with respect to
1149  /// the arc length, used in continuation
1150  inline double &dof_derivative(const unsigned &i)
1151  {
1153  {return Dof_derivative[i];}
1154  else
1155  {return (*(Dof_pt[i]+Dof_derivative_offset));}
1156 }
1157 
1158  /// \short Access function to the current value of the i-th (local)
1159  /// dof at the start of a continuation step
1160  inline double &dof_current(const unsigned &i)
1161  {
1163  {return Dof_current[i];}
1164  else
1165  {return (*(Dof_pt[i]+Dof_current_offset));}
1166  }
1167 
1168  /// \short Set initial condition (incl previous timesteps).
1169  /// We need to establish this interface
1170  /// because I.C. needs to be reset when problem is adapted during
1171  /// the first timestep.
1172  virtual void set_initial_condition()
1173  {
1174  std::ostringstream warn_message;
1175  warn_message
1176  << "Warning: We're using the default (empty) set_initial_condition().\n"
1177  << "If the initial conditions isn't re-assigned after a mesh adaption \n"
1178  << "the initial conditions will represent the interpolation of the \n"
1179  << "initial conditions that were assigned on the original coarse mesh.\n";
1180  OomphLibWarning(warn_message.str(),
1181  "Problem::set_initial_condition()",
1182  OOMPH_EXCEPTION_LOCATION);
1183 
1184  // Indicate that this function has been called. This flag is set so
1185  // that the unsteady_newton_solve routine can be instructed whether
1186  // or not to shift the history values. If set_initial_condition() has
1187  // been overloaded than this (default) function won't be called, and
1188  // so this flag will remain false (its default value). If
1189  // set_initial_condition() has not been overloaded then this function
1190  // will be called and the flag set to true.
1192  }
1193 
1194  /// \short Function to calculate a global error norm, used in adaptive
1195  /// timestepping to control the change in timestep. Individual errors for
1196  /// each data object can be obtained via the data timestepper's
1197  /// temporal_error_in_value or temporal_error_in_position functions and should
1198  /// be combined to construct a global norm. For example, in fluids problems
1199  /// a suitable norm is usually the weighted sum of the errors in the
1200  /// velocities; for moving mesh problems is it usually better to use
1201  /// the weighted sum of the errors in position.
1203  {
1204  std::string error_message =
1205  "The global_temporal_error_norm function will be problem-specific:\n";
1206  error_message +=
1207  "Please write your own in your Problem class";
1208 
1209  throw OomphLibError(error_message,
1210  OOMPH_CURRENT_FUNCTION,
1211  OOMPH_EXCEPTION_LOCATION);
1212  return 0.0;
1213  }
1214 
1215  /// The communicator for this problem
1217 
1218  public:
1219 
1220  /// access function to the oomph-lib communicator
1222  {
1223  return Communicator_pt;
1224  }
1225 
1226  /// access function to the oomph-lib communicator, const version
1228  {
1229  return Communicator_pt;
1230  }
1231 
1232  /// Function pointer for spatial error estimator
1234  Vector<double>& elemental_error);
1235 
1236  /// Function pointer for spatial error estimator with doc
1239  elemental_error,
1240  DocInfo& doc_info);
1241 
1242  /// \short Constructor: Allocate space for one time stepper
1243  /// and set all pointers to NULL and set defaults for all
1244  /// parameters.
1245  Problem();
1246 
1247  /// Broken copy constructor
1248  Problem(const Problem& dummy)
1249  {
1250  BrokenCopy::broken_copy("Problem");
1251  }
1252 
1253  /// Broken assignment operator
1254  void operator=(const Problem&)
1255  {
1256  BrokenCopy::broken_assign("Problem");
1257  }
1258 
1259 
1260  /// Virtual destructor to clean up memory
1261  virtual ~Problem();
1262 
1263  /// Return a pointer to the global mesh
1264  Mesh*& mesh_pt() {return Mesh_pt;}
1265 
1266  /// Return a pointer to the global mesh (const version)
1267  Mesh* const &mesh_pt() const {return Mesh_pt;}
1268 
1269  /// \short Return a pointer to the i-th submesh. If there are
1270  /// no submeshes, the 0-th submesh is the global mesh itself.
1271  Mesh* &mesh_pt(const unsigned &imesh)
1272  {
1273  // If there are no submeshes then the zero-th submesh is the
1274  // mesh itself
1275  if((Sub_mesh_pt.size()==0)&&(imesh==0)) {return Mesh_pt;}
1276  else {return Sub_mesh_pt[imesh];}
1277  }
1278 
1279  /// Return a pointer to the i-th submesh (const version)
1280  Mesh* const &mesh_pt(const unsigned &imesh) const
1281  {
1282  // If there are no submeshes then the zero-th submesh is the
1283  // mesh itself
1284  if((Sub_mesh_pt.size()==0)&&(imesh==0)) {return Mesh_pt;}
1285  else {return Sub_mesh_pt[imesh];}
1286  }
1287 
1288  /// Return number of submeshes
1289  unsigned nsub_mesh() const {return Sub_mesh_pt.size();}
1290 
1291  /// \short Add a submesh to the problem and return its number, i, by which it
1292  /// can be accessed via mesh_pt(i).
1293  unsigned add_sub_mesh(Mesh* const &mesh_pt)
1294  {
1295  Sub_mesh_pt.push_back(mesh_pt);
1296  return Sub_mesh_pt.size()-1;
1297  }
1298 
1299 
1300  /// \short Flush the problem's collection of sub-meshes. Must be followed
1301  /// by call to rebuild_global_mesh().
1302  void flush_sub_meshes() {Sub_mesh_pt.clear();}
1303 
1304  /// \short Build the global mesh by combining the all the submeshes.
1305  /// \b Note: The nodes boundary information refers to the boundary numbers
1306  /// within the submesh!
1307  void build_global_mesh();
1308 
1309  /// \short If one of the submeshes has changed (e.g. by
1310  /// mesh adaptation) we need to update the global mesh.
1311  /// \b Note: The nodes boundary information refers to the
1312  /// boundary numbers within the submesh!
1313  void rebuild_global_mesh();
1314 
1315 #ifdef OOMPH_HAS_MPI
1316 
1317  /// \short Function to build the Problem's base mesh; this must be supplied
1318  /// by the user if they wish to use the load_balance() functionality,
1319  /// which is only available to problems that have already been distributed.
1320  /// If the problem has multiple meshes, each mesh must be built, added as
1321  /// as a submesh, and a call to build_global_mesh() must be made
1322  /// in this function. On return from this function all meshes must
1323  /// have been refined to the same level that they were in the when
1324  /// Problem::distribute() was first called.
1325  virtual void build_mesh()
1326  {
1327  std::string error_message =
1328  "You are likely to have reached this error from Problem::load_balance()\n";
1329  error_message +=
1330  "The build_mesh() function is problem-specific and must be supplied:\n";
1331  error_message +=
1332  "Please write your own in your Problem class, ensuring that\n";
1333  error_message +=
1334  "any (sub)meshes are built before calling Problem::build_global_mesh().\n";
1335  throw OomphLibError(error_message,
1336  OOMPH_CURRENT_FUNCTION,
1337  OOMPH_EXCEPTION_LOCATION);
1338  }
1339 
1340  /// \short Balance the load of a (possibly non-uniformly refined) problem
1341  /// that has already been distributed, by re-distributing elements over
1342  /// processors.
1344  {
1345  // Dummy DocInfo
1346  DocInfo doc_info;
1347  doc_info.disable_doc();
1348 
1349  // Don't report stats
1350  bool report_stats=false;
1351 
1352  // Dummy imposed partioning vector
1353  Vector<unsigned> input_target_domain_for_local_non_halo_element;
1354 
1355  // Do it!
1356  load_balance(doc_info,report_stats,
1357  input_target_domain_for_local_non_halo_element);
1358  }
1359 
1360  /// \short Balance the load of a (possibly non-uniformly refined) problem that
1361  /// has already been distributed, by re-distributing elements over
1362  /// processors. Produce explicit stats of load balancing process if
1363  /// boolean, report_stats, is set to true.
1364  void load_balance(const bool& report_stats)
1365  {
1366  // Dummy DocInfo
1367  DocInfo doc_info;
1368  doc_info.disable_doc();
1369 
1370  // Dummy imposed partioning vector
1371  Vector<unsigned> input_target_domain_for_local_non_halo_element;
1372 
1373  // Do it
1374  load_balance(doc_info,report_stats,
1375  input_target_domain_for_local_non_halo_element);
1376  }
1377 
1378 
1379  /// \short Balance the load of a (possibly non-uniformly refined) problem that
1380  /// has already been distributed, by re-distributing elements over
1381  /// processors. Produce explicit stats of load balancing process if
1382  /// boolean, report_stats, is set to true.
1383  void load_balance(DocInfo& doc_info,
1384  const bool& report_stats)
1385  {
1386  // Dummy imposed partioning vector
1387  Vector<unsigned> input_target_domain_for_local_non_halo_element;
1388 
1389  // Do it
1390  load_balance(doc_info,report_stats,
1391  input_target_domain_for_local_non_halo_element);
1392  }
1393 
1394  /// \short Balance the load of a (possibly non-uniformly refined) problem that
1395  /// has already been distributed, by re-distributing elements over
1396  /// processors. Produce explicit stats of load balancing process if
1397  /// boolean, report_stats, is set to true and doc various bits of
1398  /// data (mainly for debugging) in directory specified by DocInfo object.
1399  /// If final input vector is non-zero-sized it provides an imposed
1400  /// partitioning.
1401  void load_balance(DocInfo& doc_info,
1402  const bool& report_stats,
1403  const Vector<unsigned>&
1404  input_target_domain_for_local_non_halo_element);
1405 
1406  /// \short Set the use of the default partition in the load balance
1409 
1410  /// \short Do not use the default partition in the load balance
1413 
1414  /// \short Load balance helper routine: refine each new base (sub)mesh
1415  /// based upon the elements to be refined within each tree at each root
1416  /// on the current processor
1418  (Vector<Vector<Vector<unsigned> > >& to_be_refined_on_each_root,
1419  const unsigned& max_level_overall);
1420 
1421 #endif
1422 
1423  /// Return a pointer to the linear solver object
1425 
1426  /// Return a pointer to the linear solver object (const version)
1428 
1429  /// Return a pointer to the linear solver object used for explicit time
1430  /// stepping.
1433 
1434  /// Return a pointer to the linear solver object used for explicit time
1435  /// stepping (const version).
1438 
1439  /// Return a pointer to the eigen solver object
1441 
1442  /// Return a pointer to the eigen solver object (const version)
1443  EigenSolver* const &eigen_solver_pt() const {return Eigen_solver_pt;}
1444 
1445  /// Return a pointer to the global time object
1446  Time* &time_pt() {return Time_pt;}
1447 
1448  /// Return a pointer to the global time object (const version).
1449  Time* time_pt() const {return Time_pt;}
1450 
1451  /// Return the current value of continuous time
1452  double& time();
1453 
1454  /// Return the current value of continuous time (const version)
1455  double time() const;
1456 
1457 
1458  /// \short Access function for the pointer to the first (presumably only)
1459  /// timestepper
1461  {
1462  if (Time_stepper_pt.size()==0)
1463  {
1464  throw OomphLibError("No timestepper allocated yet\n",
1465  OOMPH_CURRENT_FUNCTION,
1466  OOMPH_EXCEPTION_LOCATION);
1467  }
1468  return Time_stepper_pt[0];
1469  }
1470 
1471  /// \short Access function for the pointer to the first (presumably only)
1472  /// timestepper
1474  {
1475  if (Time_stepper_pt.size()==0)
1476  {
1477  throw OomphLibError("No timestepper allocated yet\n",
1478  OOMPH_CURRENT_FUNCTION,
1479  OOMPH_EXCEPTION_LOCATION);
1480  }
1481  return Time_stepper_pt[0];
1482  }
1483 
1484  /// Return a pointer to the i-th timestepper
1485  TimeStepper* &time_stepper_pt(const unsigned &i)
1486  {return Time_stepper_pt[i];}
1487 
1488  /// Return a pointer to the explicit timestepper
1490  {return Explicit_time_stepper_pt;}
1491 
1492  /// \short Set all problem data to have the same timestepper (timestepper_pt)
1493  /// Return the new number of dofs in the problem
1494  unsigned long
1496  const bool &preserve_existing_data=false);
1497 
1498  /// \short Shift all values along to prepare for next timestep
1499  virtual void shift_time_values();
1500 
1501  /// Return a pointer to the assembly handler object
1503 
1504  /// Return a pointer to the assembly handler object (const version)
1506  {return Assembly_handler_pt;}
1507 
1508  /// \short Access function to min timestep in adaptive timestepping
1509  double& minimum_dt() {return Minimum_dt;}
1510 
1511  /// \short Access function to max timestep in adaptive timestepping
1512  double& maximum_dt() {return Maximum_dt;}
1513 
1514  /// \short Access function to max Newton iterations before giving up.
1516 
1517  /// \short Access function to Problem_is_nonlinear.
1518  void problem_is_nonlinear(const bool& prob_lin)
1519  {
1520  Problem_is_nonlinear = prob_lin;
1521  }
1522 
1523 
1524  /// \short Access function to max residuals in Newton iterations before giving
1525  /// up.
1526  double& max_residuals() {return Max_residuals;}
1527 
1528  /// Access function for Time_adaptive_newton_crash_on_solve_fail.
1531 
1532  /// \short Access function to tolererance of the Newton solver, i.e. the
1533  /// maximum value of the residuals that will be accepted.
1535 
1536  /// \short Add a timestepper to the problem. The function will automatically
1537  /// create or resize the Time object so that it contains the appropriate
1538  /// number of levels of storage.
1540 
1541  /// \short Set the explicit timestepper for the problem. The function
1542  /// will automatically create or resize the Time object so that it contains
1543  /// the appropriate number of levels of storage
1545  const &explicit_time_stepper_pt);
1546 
1547 
1548  /// \short Set all timesteps to the same value, dt, and assign
1549  /// weights for all timesteppers in the problem
1550  void initialise_dt(const double& dt);
1551 
1552  /// \short Set the value of the timesteps to be equal to the values passed in
1553  /// a vector and assign weights for all timesteppers in the problem
1554  void initialise_dt(const Vector<double>& dt);
1555 
1556  /// Return a pointer to the the i-th global data object
1557  Data* &global_data_pt(const unsigned &i) {return Global_data_pt[i];}
1558 
1559  /// \short Add Data to the Problem's global data -- the Problem will
1560  /// perform equation numbering etc. for such Data
1562  {Global_data_pt.push_back(global_data_pt);}
1563 
1564 
1565  /// \short Flush the Problem's global data -- resizes container to zero.
1566  /// Data objects are not deleted!
1567  void flush_global_data() {Global_data_pt.resize(0);}
1568 
1569  /// \short Return the pointer to the dof distribution (read-only)
1571  {return Dof_distribution_pt;}
1572 
1573  /// Return the number of dofs
1574  unsigned long ndof() const
1575  {return Dof_distribution_pt->nrow();}
1576 
1577  /// Return the number of time steppers
1578  unsigned ntime_stepper() const {return Time_stepper_pt.size();}
1579 
1580  ///Return the number of global data values
1581  unsigned nglobal_data() const {return Global_data_pt.size();}
1582 
1583  /// \short Self-test: Check meshes and global data. Return 0 for OK
1584  unsigned self_test();
1585 
1586  /// \short Insist that local dof pointers are set up in each element
1587  /// when equation numbering takes place
1590 
1591  /// \short Insist that local dof pointers are NOT set up in each element
1592  /// when equation numbering takes place (the default)
1595 
1596  /// \short Assign all equation numbers for problem: Deals with global
1597  /// data (= data that isn't attached to any elements) and then
1598  /// does the equation numbering for the elements. Virtual so it
1599  /// can be overloaded in MPI problems. Bool argument can be set to false
1600  /// to ignore assigning local equation numbers (found to be necessary in
1601  /// the parallel implementation of locate_zeta between multiple meshes).
1602  unsigned long assign_eqn_numbers(const bool&
1603  assign_local_eqn_numbers=true);
1604 
1605  /// \short Function to describe the dofs in terms of the global
1606  /// equation number, i.e. what type of value (nodal value of
1607  /// a Node; value in a Data object; value of internal Data in an
1608  /// element; etc) is the unknown with a certain global equation number.
1609  /// Output stream defaults to oomph_info.
1610  void describe_dofs(std::ostream& out=*(oomph_info.stream_pt())) const;
1611 
1612  /// \short Indicate that the problem involves discontinuous elements
1613  /// This allows for a more efficiently assembly and inversion of the
1614  /// mass matrix
1617 
1618  /// \short Disable the use of a discontinuous-element formulation.
1619  /// Note that the methods will all still work even if the elements are
1620  /// discontinuous, we will just be assembling a larger system matrix than
1621  /// necessary.
1624 
1625  /// \short Return the vector of dofs, i.e. a vector containing the current
1626  /// values of all unknowns.
1627  void get_dofs(DoubleVector& dofs) const;
1628 
1629  /// \short Return vector of the t'th history value of all dofs.
1630  void get_dofs(const unsigned& t, DoubleVector& dofs) const;
1631 
1632  /// \short Set the values of the dofs
1633  void set_dofs(const DoubleVector &dofs);
1634 
1635  /// \short Set the history values of the dofs
1636  void set_dofs(const unsigned& t, DoubleVector& dofs);
1637 
1638  /// Set history values of dofs from the type of vector stored in
1639  /// problem::Dof_pt.
1640  void set_dofs(const unsigned& t, Vector<double*>& dof_pt);
1641 
1642  /// \short Add lambda x incremenet_dofs[l] to the l-th dof
1643  void add_to_dofs(const double &lambda, const DoubleVector &increment_dofs);
1644 
1645  /// \short Return a pointer to the dof, indexed by global equation number
1646  /// which may be haloed or stored locally. If it is haloed, a local copy
1647  /// must have been requested on this processor in the Halo_scheme_pt.
1648  inline double* global_dof_pt(const unsigned &i)
1649  {
1650 #ifdef OOMPH_HAS_MPI
1651  //If the problem is distributed I have to do something different
1653  {
1654 #ifdef PARANOID
1655  if(Halo_scheme_pt==0)
1656  {
1657  std::ostringstream error_stream;
1658  error_stream
1659  << "Direct access to global dofs in a distributed problem is only\n"
1660  << "possible if the function setup_dof_halo_scheme() has been called.\n"
1661  << "You can access local dofs via Problem::dof(i).\n";
1662  throw OomphLibError(error_stream.str(),
1663  OOMPH_CURRENT_FUNCTION,
1664  OOMPH_EXCEPTION_LOCATION);
1665  }
1666 #endif
1667 
1668  //Work out the local coordinate of the dof
1669  //based on the original distribution stored in the Halo_scheme
1670  const unsigned first_row_local =
1672  const unsigned n_row_local =
1674 
1675  //If we are in range then just call the local value
1676  if((i >= first_row_local) && (i < first_row_local + n_row_local))
1677  {
1678  return this->Dof_pt[i-first_row_local];
1679  }
1680  //Otherwise the entry is not stored in the local processor
1681  //and we must have haloed it
1682  else
1683  {
1685  }
1686  }
1687  //Otherwise just return the dof
1688  else
1689 #endif
1690  {
1691  return this->Dof_pt[i];
1692  }
1693  }
1694 
1695  /// \short i-th dof in the problem
1696  double& dof(const unsigned& i)
1697  {
1698  return *(Dof_pt[i]);
1699  }
1700 
1701  /// \short i-th dof in the problem (const version)
1702  double dof(const unsigned& i) const
1703  {
1704  return *(Dof_pt[i]);
1705  }
1706 
1707  /// \short Pointer to i-th dof in the problem
1708  double* &dof_pt(const unsigned &i) {return Dof_pt[i];}
1709 
1710 /// \short Pointer to i-th dof in the problem (const version)
1711  double* dof_pt(const unsigned &i) const {return Dof_pt[i];}
1712 
1713  ///\short Return the residual vector multiplied by the inverse mass matrix
1714  ///Virtual so that it can be overloaded for mpi problems
1716 
1717  /// \short Get the time derivative of all values (using
1718  /// get_inverse_mass_matrix_times_residuals(..) with all time steppers set
1719  /// to steady) e.g. for use in explicit time steps. The approach used is
1720  /// slighty hacky, beware if you have a residual which is non-linear or
1721  /// implicit in the derivative or if you have overloaded
1722  /// get_jacobian(...).
1723  virtual void get_dvaluesdt(DoubleVector &f);
1724 
1725  /// \short Return the fully-assembled residuals Vector for the problem:
1726  /// Virtual so it can be overloaded in for mpi problems
1727  virtual void get_residuals(DoubleVector &residuals);
1728 
1729  /// \short Return the fully-assembled Jacobian and residuals for the problem
1730  /// Interface for the case when the Jacobian matrix is dense.
1731  /// This is virtual so, if we feel like it (e.g. for testing iterative
1732  /// solvers with specific test matrices, we can bypass the default
1733  /// assembly procedure for the Jacobian and the residual vector.
1734  virtual void get_jacobian(DoubleVector &residuals,
1735  DenseDoubleMatrix& jacobian);
1736 
1737  /// \short Return the fully-assembled Jacobian and residuals for the problem.
1738  /// Interface for the case when the Jacobian is in row-compressed storage
1739  /// format. This is virtual so, if we feel like it (e.g. for testing iterative
1740  /// solvers with specific test matrices), we can bypass the default
1741  /// assembly procedure for the Jacobian and the residual vector.
1742  virtual void get_jacobian(DoubleVector &residuals,
1743  CRDoubleMatrix &jacobian);
1744 
1745  /// \short Return the fully-assembled Jacobian and residuals for the problem.
1746  /// Interface for the case when the Jacobian is in column-compressed storage
1747  /// format. This is virtual so, if we feel like it (e.g. for testing iterative
1748  /// solvers with specific test matrices), we can bypass the default
1749  /// assembly procedure for the Jacobian and the residual vector.
1750  virtual void get_jacobian(DoubleVector &residuals,
1751  CCDoubleMatrix &jacobian);
1752 
1753  /// \short Dummy virtual function that must be overloaded by the problem to
1754  /// specify which matrices should be summed to give the final Jacobian.
1755  virtual void get_jacobian(DoubleVector &residuals,
1756  SumOfMatrices &jacobian)
1757  {
1758  std::ostringstream error_stream;
1759  error_stream
1760  << "You must overload this function in your problem class to specify\n"
1761  << "which matrices should be summed to give the final Jacobian.";
1762  throw OomphLibError(error_stream.str(),
1763  OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1764  }
1765 
1766  /// \short Return the fully-assembled Jacobian and residuals, generated by
1767  /// finite differences
1768  void get_fd_jacobian(DoubleVector &residuals,
1769  DenseMatrix<double> &jacobian);
1770 
1771  /// \short Get the derivative of the entire residuals vector wrt a
1772  /// global parameter, used in continuation problems
1773  void get_derivative_wrt_global_parameter(double* const &parameter_pt,
1774  DoubleVector &result);
1775 
1776  /// \short Return the product of the global hessian (derivative of Jacobian
1777  /// matrix with respect to all variables) with
1778  /// an eigenvector, Y, and any number of other specified vectors C
1779  /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}.
1780  /// This function is used in assembling and solving the augmented systems
1781  /// associated with bifurcation tracking.
1782  /// The default implementation is to use finite differences at the global
1783  /// level.
1785  DoubleVectorWithHaloEntries const &Y,
1788 
1789 
1790 
1791  /// \short Get derivative of an element in the problem wrt a global parameter,
1792  /// used in continuation problems
1793  //void get_derivative_wrt_global_parameter(double* const &parameter_pt,
1794  // GeneralisedElement* const &elem_pt,
1795  // Vector<double> &result);
1796 
1797  /// \short Solve an eigenproblem as assembled by EigenElements
1798  /// calculate n_eval eigenvalues and return the corresponding
1799  /// eigenvectors. The boolean flag (default true) specifies whether
1800  /// the steady jacobian should be assembled. If the flag is false
1801  /// then the weighted mass-matrix terms from the timestepper will
1802  /// be included in the jacobian --- this is almost certainly never
1803  /// wanted.
1804  void solve_eigenproblem(const unsigned &n_eval,
1805  Vector<std::complex<double> > &eigenvalue,
1806  Vector<DoubleVector> &eigenvector,
1807  const bool &steady=true);
1808 
1809  /// \short Solve an eigenproblem as assembled by EigenElements,
1810  /// but only return the eigenvalues, not the eigenvectors.
1811  /// The boolean flag (default true) is used to specify whether the
1812  /// weighted mass-matrix terms from the timestepping scheme should
1813  /// be included in the jacobian.
1814  void solve_eigenproblem(const unsigned &n_eval,
1815  Vector<std::complex<double> > &eigenvalue,
1816  const bool &steady=true)
1817  {
1818  //Create temporary storage for the eigenvectors (potentially wasteful)
1819  Vector<DoubleVector> eigenvector;
1820  solve_eigenproblem(n_eval,eigenvalue,eigenvector,steady);
1821  }
1822 
1823  /// \short Get the matrices required by a eigensolver. If the
1824  ///shift parameter is non-zero the second matrix will be shifted
1825  virtual void get_eigenproblem_matrices(CRDoubleMatrix &mass_matrix,
1826  CRDoubleMatrix &main_matrix,
1827  const double &shift=0.0);
1828 
1829  /// \short Assign the eigenvector passed to the function
1830  /// to the dofs in the problem so that it can be output by
1831  /// the usual routines.
1832  void assign_eigenvector_to_dofs(DoubleVector &eigenvector);
1833 
1834  /// \short Add the eigenvector passed to the function scaled by
1835  /// the constat epsilon to the dofs in the problem so that
1836  /// it can be output by the usual routines
1837  void add_eigenvector_to_dofs(const double &epsilon,
1838  const DoubleVector &eigenvector);
1839 
1840 
1841  /// \short Store the current values of the degrees of freedom
1842  void store_current_dof_values();
1843 
1844  /// \short Restore the stored values of the degrees of freedom
1845  void restore_dof_values();
1846 
1847  /// \short Enable recycling of Jacobian in Newton iteration
1848  /// (if the linear solver allows it).
1849  /// Useful for linear problems with constant Jacobians or nonlinear
1850  /// problems where you're willing to risk the trade-off between
1851  /// faster solve times and degraded Newton convergence rate
1853  {
1856  }
1857 
1858  /// \short Disable recycling of Jacobian in Newton iteration
1860  {
1863  }
1864 
1865  /// \short Is recycling of Jacobian in Newton iteration enabled?
1867  {
1869  }
1870 
1872  {
1874  }
1875 
1876  /// \short Use Newton method to solve the problem
1877  void newton_solve();
1878 
1879  /// \short enable globally convergent Newton method
1881  {
1883  }
1884 
1885  /// \short disable globally convergent Newton method
1887  {
1889  }
1890 
1891  private:
1892 
1893  /// Line search helper for globally convergent Newton method
1895  const double& half_residual_squared_old,
1896  DoubleVector& gradient,
1897  DoubleVector& newton_dir,
1898  double& half_residual_squared,
1899  const double& stpmax);
1900 
1901  public:
1902 
1903  /// \short Adaptive Newton solve: up to max_adapt adaptations of all
1904  /// refineable submeshes are performed to achieve the
1905  /// the error targets specified in the refineable submeshes.
1906  void newton_solve(unsigned const &max_adapt);
1907 
1908  /// \short Solve a steady problem using adaptive Newton's method,
1909  /// but in the context of an overall unstady problem, perhaps to
1910  /// determine an initial condition.
1911  /// This is achieved by setting the weights in the timesteppers to be zero
1912  /// which has the effect of rendering them steady timesteppers.
1913  /// The optional argument max_adapt specifies the max. number of
1914  /// adaptations of all refineable submeshes are performed to
1915  /// achieve the the error targets specified in the refineable submeshes.
1916  void steady_newton_solve(unsigned const &max_adapt=0);
1917 
1918 /// \short Copy Data values, nodal positions etc from specified problem.
1919 /// Note: This is not a copy constructor. We assume that the current
1920 /// and the "original" problem have both been created by calling
1921 /// the same problem constructor so that all Data objects,
1922 /// time steppers etc. in the two problems are completely independent.
1923 /// This function copies the nodal, internal and global values,
1924 /// and the time parameters from the original problem into
1925 /// "this" one. This functionality is required, e.g. for
1926 /// multigrid computations.
1927  void copy(Problem* orig_problem_pt);
1928 
1929  /// \short Make and return a pointer to the copy of the problem. A virtual
1930  /// function that must be filled in by the user is they wish to perform
1931  /// adaptive refinement in bifurcation tracking or in multigrid problems.
1932  /// ALH: WILL NOT BE NECESSARY IN BIFURCATION TRACKING IN LONG RUN...
1933  virtual Problem* make_copy();
1934 
1935  /// \short Read refinement pattern of all refineable meshes and refine them
1936  /// accordingly, then read all Data and nodal position info from
1937  /// file for restart. Return flag to indicate if the restart was from
1938  /// steady or unsteady solution
1939  virtual void read(std::ifstream& restart_file, bool& unsteady_restart);
1940 
1941  /// \short Read refinement pattern of all refineable meshes and refine them
1942  /// accordingly, then read all Data and nodal position info from
1943  /// file for restart.
1944  virtual void read(std::ifstream& restart_file)
1945  {
1946  bool unsteady_restart;
1947  read(restart_file,unsteady_restart);
1948  }
1949 
1950  /// \short Dump refinement pattern of all refineable meshes and all generic
1951  /// Problem data to file for restart.
1952  virtual void dump(std::ofstream& dump_file) const;
1953 
1954  /// \short Dump refinement pattern of all refineable meshes and all generic
1955  /// Problem data to file for restart.
1956  void dump(const std::string &dump_file_name) const
1957  {
1958  std::ofstream dump_stream(dump_file_name.c_str());
1959 #ifdef PARANOID
1960  if(!dump_stream.is_open())
1961  {
1962  std::string err = "Couldn't open file "+dump_file_name;
1963  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
1964  OOMPH_CURRENT_FUNCTION);
1965  }
1966 #endif
1967  dump(dump_stream);
1968  }
1969 
1970 #ifdef OOMPH_HAS_MPI
1971 
1972  /// \short Get pointers to all possible halo data indexed by global
1973  /// equation number in a map.
1974  void get_all_halo_data(std::map<unsigned,double*> &map_of_halo_data);
1975 
1976  /// \short Classify any non-classified nodes into halo/haloed and
1977  /// synchronise equation numbers. Return the total
1978  /// number of degrees of freedom in the overall problem
1979  long synchronise_eqn_numbers(const bool& assign_local_eqn_numbers=true);
1980 
1981  /// \short Synchronise the degrees of freedom by overwriting
1982  /// the haloed values with their non-halo counterparts held
1983  /// on other processors. Bools control if we deal with data associated with
1984  /// external halo/ed elements/nodes or the "normal" halo/ed ones.
1985  void synchronise_dofs(const bool& do_halos,
1986  const bool& do_external_halos);
1987 
1988  /// \short Perform all required synchronisation in solvers
1989  void synchronise_all_dofs();
1990 
1991  /// Check the halo/haloed node/element schemes
1992  void check_halo_schemes(DocInfo& doc_info);
1993 
1994  /// Check the halo/haloed node/element schemes
1996  {
1997  DocInfo tmp_doc_info;
1998  tmp_doc_info.disable_doc();
1999  check_halo_schemes(tmp_doc_info);
2000  }
2001 
2002  /// \short Distribute the problem and doc, using the specified partition;
2003  /// returns a vector which details the partitioning
2004  Vector<unsigned> distribute(const Vector<unsigned>& element_partition,
2005  DocInfo& doc_info,
2006  const bool& report_stats=false);
2007 
2008  /// \short Distribute the problem; returns a vector which
2009  /// details the partitioning
2011  const bool& report_stats=false);
2012 
2013  /// \short Distribute the problem using the specified partition;
2014  /// returns a vector which details the partitioning
2015  Vector<unsigned> distribute(const Vector<unsigned>& element_partition,
2016  const bool& report_stats=false);
2017 
2018  /// \short Distribute the problem; returns a vector which
2019  /// details the partitioning
2020  Vector<unsigned> distribute(const bool& report_stats=false);
2021 
2022  /// /short Partition the global mesh, return vector specifying the processor
2023  /// number for each element. Virtual so that it can be overloaded by
2024  /// any user; the default is to use METIS to perform the partitioning
2025  /// (with a bit of cleaning up afterwards to sort out "special cases").
2026  virtual void partition_global_mesh(Mesh* &global_mesh_pt, DocInfo& doc_info,
2027  Vector<unsigned>& element_domain,
2028  const bool& report_stats=false);
2029 
2030  /// \short (Irreversibly) prune halo(ed) elements and nodes, usually
2031  /// after another round of refinement, to get rid of
2032  /// excessively wide halo layers. Note that the current
2033  /// mesh will be now regarded as the base mesh and no unrefinement
2034  /// relative to it will be possible once this function
2035  /// has been called.
2036  void prune_halo_elements_and_nodes(DocInfo& doc_info,
2037  const bool& report_stats);
2038 
2039  /// \short (Irreversibly) prune halo(ed) elements and nodes, usually
2040  /// after another round of refinement, to get rid of
2041  /// excessively wide halo layers. Note that the current
2042  /// mesh will be now regarded as the base mesh and no unrefinement
2043  /// relative to it will be possible once this function
2044  /// has been called.
2045  void prune_halo_elements_and_nodes(const bool& report_stats=false)
2046  {
2047  DocInfo doc_info;
2048  doc_info.disable_doc();
2049  prune_halo_elements_and_nodes(doc_info,report_stats);
2050  }
2051 
2052  /// Threshold for error throwing in Problem::check_halo_schemes()
2054 
2055  /// Access to Problem_has_been_distributed flag
2057  {
2059  }
2060 
2061 #endif
2062 
2063  /// \short Wrapper function to delete external storage for
2064  /// each submesh of the problem
2066 
2067  /// \short Boolean to indicate if all output is suppressed in
2068  /// Problem::newton_solve(). Defaults to false.
2070 
2071 
2072  protected:
2073 
2074  /// \short Boolean to indicate whether a Newton step should be taken
2075  /// even if the initial residuals are below the required tolerance
2077 
2078  /// \short What it says: If temporally adaptive Newton solver fails to
2079  /// to converge, reduce timestep by this factor and try again; defaults
2080  /// to 1/2; can be over-written by user in derived problem.
2082 
2083  /// \short Perform a basic arc-length continuation step using Newton's
2084  /// method. Returns number of Newton steps taken.
2085  unsigned newton_solve_continuation(double* const &parameter_pt);
2086 
2087  /// \short This function performs a basic continuation step using the Newton
2088  /// method. The number of Newton steps taken is returned, to be used in any
2089  /// external step-size control routines.
2090  /// The governing parameter of the problem is passed as a pointer to the
2091  /// routine, as is a vector in which
2092  /// to store the derivatives wrt the parameter, if required.
2093  unsigned newton_solve_continuation(double* const &parameter_pt,
2094  DoubleVector &z);
2095 
2096  ///\short A function to calculate the derivatives wrt the arc-length. This
2097  /// version of the function actually does a linear solve so that the
2098  /// derivatives
2099  /// are calculated "exactly" rather than using the values at the Newton
2100  /// step just before convergence. This is necessary in spatially adaptive
2101  /// problems, in which the number of degrees of freedom changes and so
2102  /// the appropriate derivatives must be calculated for the new variables.
2103  /// This function is called if no Newton steps were taken in the
2104  /// continuation routine ... i.e. the initial residuals were sufficiently
2105  /// small.
2106  void calculate_continuation_derivatives(double* const &parameter_pt);
2107 
2108  /// \short A function to calculate the derivatives with respect to the
2109  /// arc-length required for continuation. The arguments is the solution
2110  /// of the linear system,
2111  /// Jz = dR/dparameter, that gives du/dparameter and the direction
2112  /// output from the newton_solve_continuation function. The derivatives
2113  /// are stored in the ContinuationParameters namespace.
2115 
2116  /// \short A function to calculate the derivatives with respect to the
2117  /// arc-length required for continuation by finite differences, using
2118  /// the previous values of the solution. The derivatives are stored in
2119  /// the ContinuationParameters namespace.
2120  void calculate_continuation_derivatives_fd(double* const &parameter_pt);
2121 
2122  /// \short Return a boolean flag to indicate whether the pointer
2123  /// parameter_pt refers to values stored in a Data object that
2124  /// is contained within the problem
2125  bool does_pointer_correspond_to_problem_data(double* const &parameter_pt);
2126 
2127  public:
2128 
2129  /// \short Virtual function that is used to symmetrise the problem so that
2130  /// the current solution exactly satisfies any symmetries within the system.
2131  /// Used when adpativly solving pitchfork detection problems when small
2132  /// asymmetries in the coarse solution can be magnified
2133  /// leading to very inaccurate answers on the fine mesh.
2134  /// This is always problem-specific and must be filled in by the user
2135  /// The default issues a warning
2137 
2138  ///\short Return pointer to the parameter that is used in the
2139  /// bifurcation detection. If we are not tracking a bifurcation then
2140  /// an error will be thrown by the AssemblyHandler
2141  double* bifurcation_parameter_pt() const;
2142 
2143 
2144  /// \short Return the eigenfunction calculated as part of a
2145  /// bifurcation tracking process. If we are not tracking a bifurcation
2146  /// then an error will be thrown by the AssemblyHandler
2148 
2149 
2150  /// \short Turn on fold tracking using the augmented system specified
2151  /// in the FoldHandler class. After a call to this function subsequent calls
2152  /// of the standard solution methods will converge to a fold (limit) point
2153  /// at a particular value of the variable addressed by parameter_pt.
2154  /// The system may not converge if the initial guess is sufficiently poor
2155  /// or, alternatively, if finite differencing is used to calculate the
2156  /// jacobian matrix in the elements. If the boolean flag block_solver is
2157  /// true (the default) then a block factorisation is used to solve the
2158  /// augmented system which is both faster and uses less memory.
2159  void activate_fold_tracking(double* const &parameter_pt,
2160  const bool &block_solve=true);
2161 
2162  ///\short Activate generic bifurcation tracking for a single (real) eigenvalue
2163  ///where the initial guess for the eigenvector can be specified.
2165  double* const &parameter_pt,
2166  const DoubleVector &eigenvector,
2167  const bool &block_solve=true);
2168 
2169  ///\short Activate generic bifurcation tracking for a single (real) eigenvalue
2170  ///where the initial guess for the eigenvector can be specified and the
2171  ///normalisation condition can also be specified.
2173  double* const &parameter_pt,
2174  const DoubleVector &eigenvector,
2175  const DoubleVector &normalisation,
2176  const bool &block_solve=true);
2177 
2178 
2179  /// \short Turn on pitchfork tracking using the augmented system specified
2180  /// in the PitchForkHandler class.
2181  /// After a call to this function subsequent calls
2182  /// of the standard solution methods will converge to a pitchfork
2183  /// bifurcation
2184  /// at a particular value of the variable addressed by parameter_pt.
2185  /// The symmetry that is to be broken must be specified by supplying
2186  /// a symmetry_vector(ndof). The easiest way to determine such a vector
2187  /// is to solve the associated eigenproblem \f$ Jx = \lambda M x\f$
2188  /// and pass in the eigenvector. This is not always necessary however,
2189  /// if the symmetry is easy to construct.
2190  /// The system may not converge if the initial guess is sufficiently poor
2191  /// or, alternatively, if finite differencing is used to calculate the
2192  /// jacobian matrix in the elements. If the boolean flag block_solver is
2193  /// true (the default) then a block factorisation is used to solve the
2194  /// augmented system which is both faster and requires less memory.
2195  void activate_pitchfork_tracking(double* const &parameter_pt,
2196  const DoubleVector &symmetry_vector,
2197  const bool &block_solve=true);
2198 
2199  /// \short Turn on Hopf bifurcation
2200  /// tracking using the augmented system specified
2201  /// in the HopfHandler class. After a call to this function subsequent calls
2202  /// of the standard solution methods will converge to a Hopf bifuraction
2203  /// at a particular value of the variable addressed by parameter_pt.
2204  /// The system may not converge if the initial guess is sufficiently poor
2205  /// or, alternatively, if finite differencing is used to calculate the
2206  /// jacobian matrix in the elements.
2207  void activate_hopf_tracking(double* const &parameter_pt,
2208  const bool &block_solve=true);
2209 
2210  /// \short Turn on Hopf bifurcation
2211  /// tracking using the augmented system specified
2212  /// in the HopfHandler class. After a call to this function subsequent calls
2213  /// of the standard solution methods will converge to a Hopf bifuraction
2214  /// at a particular value of the variable addressed by parameter_pt.
2215  /// The system may not converge if the initial guess is sufficiently poor
2216  /// or, alternatively, if finite differencing is used to calculate the
2217  /// jacobian matrix in the elements. This interface allows specification
2218  /// of an inital guess for the frequency and real and imaginary parts
2219  /// of the null vector, such as might be obtained from an eigensolve
2220  void activate_hopf_tracking(double* const &parameter_pt,
2221  const double &omega,
2222  const DoubleVector &null_real,
2223  const DoubleVector &null_imag,
2224  const bool &block_solve=true);
2225 
2226  /// \short Deactivate all bifuraction tracking, by reseting
2227  /// the assembly handler to the default
2230 
2231  /// \short Reset the system to the standard non-augemented state
2233 
2234  private:
2235 
2236  /// \short Private helper function that actually contains the guts
2237  /// of the arc-length stepping, parameter_pt is a pointer to the
2238  /// parameter that is traded for the arc-length constraint, ds is
2239  /// the desired arc length and max_adapt is the maximum number of
2240  /// spatial adaptations. The pointer to the parameter may be changed
2241  /// if this is called from the Data-based interface
2242  double arc_length_step_solve_helper(double* const &parameter_pt,
2243  const double &ds,
2244  const unsigned &max_adapt);
2245 
2246 
2247  //ALH_DEVELOP
2248  protected:
2249 
2250  /// \short Private helper function that is used to set the appropriate
2251  /// pinned values for continuation.
2253 
2254  public:
2255 
2256  /// \short Solve a steady problem using arc-length continuation, when the
2257  /// parameter that becomes a variable corresponding to the arc-length
2258  /// constraint equation is an external double:
2259  /// parameter_pt is a pointer to that double,
2260  /// ds is the desired arc_length and max_adapt is the maximum
2261  /// number of spatial adaptations (default zero, no adaptation).
2262  double arc_length_step_solve(double* const &parameter_pt,
2263  const double &ds,
2264  const unsigned &max_adapt=0);
2265 
2266  /// \short Solve a steady problem using arc-length continuation, when the
2267  /// variable corresponding to the arc-length
2268  /// constraint equation is already stored in data used in the problem:
2269  /// data_pt is a pointer to the appropriate data object,
2270  /// data_index is the index of the value that will be traded for the
2271  /// constriant,
2272  /// ds is the desired arc_length and max_adapt is the maximum
2273  /// number of spatial adaptations (default zero, no adaptation).
2274  /// Note that the value must be pinned in order for this formulation to work
2275  double arc_length_step_solve(Data* const &data_pt,
2276  const unsigned &data_index,
2277  const double &ds,
2278  const unsigned &max_adapt=0);
2279 
2280 
2281 
2282  /// \short Reset the "internal" arc-length continuation parameters, so as
2283  /// to allow continuation in another parameter. N.B. The parameters that
2284  /// are reset are the "minimum" that are required, others should perhaps
2285  /// be reset, depending upon the application.
2287  {
2288  Theta_squared = 1.0;
2289  Sign_of_jacobian = 0;
2290  Continuation_direction = 1.0;
2291  Parameter_derivative = 1.0;
2293  Arc_length_step_taken=false;
2294  Dof_derivative.resize(0);
2295  }
2296 
2297  /// \short Access function for the sign of the global jacobian matrix.
2298  /// This will be set by the linear solver, if possible (direct solver).
2299  /// If not alternative methods must be used to detect bifurcations
2300  /// (solving the associated eigenproblem).
2302 
2303  /// \short Take an explicit timestep of size dt and optionally shift
2304  /// any stored values of the time history
2305  void explicit_timestep(const double &dt, const bool &shift_values=true);
2306 
2307  /// \short Advance time by dt and solve by Newton's method.
2308  /// This version always shifts time values
2309  void unsteady_newton_solve(const double &dt);
2310 
2311  /// \short Advance time by dt and solve the system,
2312  /// using Newton's method. The boolean flag is used to control
2313  /// whether the time
2314  /// values should be shifted. If it is true the current data values will
2315  /// be shifted (copied to the locations where there are stored as
2316  /// previous timesteps) before solution.
2317  void unsteady_newton_solve(const double &dt, const bool &shift_values);
2318 
2319  /// \short Unsteady adaptive Newton solve: up to max_adapt adaptations of all
2320  /// refineable submeshes are performed to achieve the
2321  /// the error targets specified in the refineable submeshes.
2322  /// If first==true, the initial conditions
2323  /// are re-assigned after the mesh adaptations.
2324  /// Shifting of time can be suppressed by overwriting the
2325  /// default value of shift (true). [Shifting must be done
2326  /// if first_timestep==true because we're constantly re-assigning
2327  /// the initial conditions; if first_timestep==true and shift==false
2328  /// shifting is performed anyway and a warning is issued.
2329  void unsteady_newton_solve(const double &dt,
2330  const unsigned &max_adapt,
2331  const bool &first,
2332  const bool& shift=true);
2333 
2334 
2335  /// \short Unsteady "doubly" adaptive Newton solve: Does temporal
2336  /// adaptation first, i.e. we try to do a timestep with an increment
2337  /// of dt, and adjusting dt until the solution on the given mesh satisfies
2338  /// the temporal error measure with tolerance epsilon. Following
2339  /// this, we do up to max_adapt spatial adaptions (without
2340  /// re-examining the temporal error). If first==true, the initial conditions
2341  /// are re-assigned after the mesh adaptations.
2342  /// Shifting of time can be suppressed by overwriting the
2343  /// default value of shift (true). [Shifting must be done
2344  /// if first_timestep==true because we're constantly re-assigning
2345  /// the initial conditions; if first_timestep==true and shift==false
2346  /// shifting is performed anyway and a warning is issued.
2347  double doubly_adaptive_unsteady_newton_solve(const double &dt,
2348  const double &epsilon,
2349  const unsigned &max_adapt,
2350  const bool &first,
2351  const bool& shift=true)
2352  {
2353  // Call helper function with default setting (do re-solve after
2354  // spatial adaptation)
2355  unsigned suppress_resolve_after_spatial_adapt_flag=0;
2357  (dt,
2358  epsilon,
2359  max_adapt,
2360  suppress_resolve_after_spatial_adapt_flag,
2361  first,
2362  shift);
2363  }
2364 
2365 
2366  /// \short Unsteady "doubly" adaptive Newton solve: Does temporal
2367  /// adaptation first, i.e. we try to do a timestep with an increment
2368  /// of dt, and adjusting dt until the solution on the given mesh satisfies
2369  /// the temporal error measure with tolerance epsilon. Following
2370  /// this, we do up to max_adapt spatial adaptions (without
2371  /// re-examining the temporal error). If first==true, the initial conditions
2372  /// are re-assigned after the mesh adaptations.
2373  /// Shifting of time can be suppressed by overwriting the
2374  /// default value of shift (true). [Shifting must be done
2375  /// if first_timestep==true because we're constantly re-assigning
2376  /// the initial conditions; if first_timestep==true and shift==false
2377  /// shifting is performed anyway and a warning is issued.
2378  /// Pseudo-Boolean flag suppress_resolve_after_spatial_adapt [0: false;
2379  /// 1: true] does what it says.]
2381  const double &dt,
2382  const double &epsilon,
2383  const unsigned &max_adapt,
2384  const unsigned& suppress_resolve_after_spatial_adapt_flag,
2385  const bool &first,
2386  const bool& shift=true)
2387  {
2388  // Call helper function
2390  (dt,
2391  epsilon,
2392  max_adapt,
2393  suppress_resolve_after_spatial_adapt_flag,
2394  first,
2395  shift);
2396  }
2397 
2398 
2399  /// \short Attempt to advance timestep by dt_desired. If the solution fails
2400  /// the timestep will be halved until convergence is achieved, or the timestep
2401  /// falls below NewtonSolverParameters::Minimum_time_step. The error control
2402  /// parameter epsilon represents the (approximate) desired magnitude of the
2403  /// global error at each timestep. The routine returns
2404  /// a double that is the suggested next timestep and should be passed as
2405  /// dt_desired the next time the routine is called. This version always
2406  /// shifts the time values.
2407  double adaptive_unsteady_newton_solve(const double &dt_desired,
2408  const double &epsilon);
2409 
2410  /// \short Attempt to advance timestep by dt_desired. If the solution fails
2411  /// the timestep will be halved until convergence is achieved, or the timestep
2412  /// falls below NewtonSolverParameters::Minimum_time_step. The error control
2413  /// parameter epsilon represents the (approximate) desired magnitude of the
2414  /// global error at each timestep. The routine returns
2415  /// a double that is the suggested next timestep and should be passed as
2416  /// dt_desired the next time the routine is called.
2417  /// Once again the boolean flag, shift_values, is used to control whether
2418  /// the time values are shifted.
2419  double adaptive_unsteady_newton_solve(const double &dt_desired,
2420  const double &epsilon,
2421  const bool &shift_values);
2422 
2423  /// \short Initialise data and nodal positions to simulate impulsive
2424  /// start from initial configuration/solution
2426 
2427  /// \short Initialise data and nodal positions to simulate an impulsive
2428  /// start and also set the initial and previous values of dt
2429  void assign_initial_values_impulsive(const double &dt);
2430 
2431  /// \short Calculate predictions
2432  void calculate_predictions();
2433 
2434  ///\short Enable recycling of the mass matrix in explicit timestepping
2435  ///schemes. Useful for timestepping on fixed meshes when you want
2436  ///to avoid the linear solve phase.
2437  void enable_mass_matrix_reuse();
2438 
2439  ///\short Turn off recyling of the mass matrix in explicit timestepping
2440  ///schemes
2442 
2443  ///\short Return whether the mass matrix is being reused
2445 
2446  /// \short Refine refineable sub-meshes, each as many times as
2447  /// specified in the vector and rebuild problem
2448  void refine_uniformly(const Vector<unsigned>& nrefine_for_mesh)
2449  {
2450  DocInfo doc_info;
2451  doc_info.disable_doc();
2452  bool prune=false;
2453  refine_uniformly_aux(nrefine_for_mesh,doc_info,prune);
2454  }
2455 
2456  /// \short Refine refineable sub-meshes, each as many times as
2457  /// specified in the vector and rebuild problem; doc refinement process
2458  void refine_uniformly(const Vector<unsigned>& nrefine_for_mesh,
2459  DocInfo& doc_info)
2460  {
2461  bool prune=false;
2462  refine_uniformly_aux(nrefine_for_mesh,doc_info,prune);
2463  }
2464 
2465  /// \short Refine refineable sub-meshes, each as many times as
2466  /// specified in the vector and rebuild problem. Prune after
2467  /// refinements
2468  void refine_uniformly_and_prune(const Vector<unsigned>& nrefine_for_mesh)
2469  {
2470  DocInfo doc_info;
2471  doc_info.disable_doc();
2472  bool prune=true;
2473  refine_uniformly_aux(nrefine_for_mesh,doc_info,prune);
2474  }
2475 
2476  /// \short Refine refineable sub-meshes, each as many times as
2477  /// specified in the vector and rebuild problem; doc refinement process
2478  void refine_uniformly_and_prune(const Vector<unsigned>& nrefine_for_mesh,
2479  DocInfo& doc_info)
2480  {
2481  bool prune=true;
2482  refine_uniformly_aux(nrefine_for_mesh,doc_info,prune);
2483  }
2484 
2485  /// \short Refine (all) refineable (sub)mesh(es) uniformly and
2486  /// rebuild problem; doc refinement process.
2487  void refine_uniformly(DocInfo& doc_info)
2488  {
2489  // Number of (sub)meshes
2490  unsigned nmesh=std::max(unsigned(1),nsub_mesh());
2491 
2492  // Refine each mesh once
2493  Vector<unsigned> nrefine_for_mesh(nmesh,1);
2494  refine_uniformly(nrefine_for_mesh);
2495  }
2496 
2497 
2498  /// \short Refine (all) refineable (sub)mesh(es) uniformly and
2499  /// rebuild problem; doc refinement process.
2501  {
2502  // Number of (sub)meshes
2503  unsigned nmesh=std::max(unsigned(1),nsub_mesh());
2504 
2505  // Refine each mesh once
2506  Vector<unsigned> nrefine_for_mesh(nmesh,1);
2507  refine_uniformly_and_prune(nrefine_for_mesh);
2508  }
2509 
2510 
2511  /// \short Refine (all) refineable (sub)mesh(es) uniformly and
2512  /// rebuild problem
2514  {
2515  DocInfo doc_info;
2516  doc_info.disable_doc();
2517  refine_uniformly(doc_info);
2518  }
2519 
2520  /// Do uniform refinement for submesh i_mesh with documentation
2521  void refine_uniformly(const unsigned& i_mesh, DocInfo& doc_info);
2522 
2523  /// Do uniform refinement for submesh i_mesh without documentation
2524  void refine_uniformly(const unsigned& i_mesh)
2525  {
2526  DocInfo doc_info;
2527  doc_info.disable_doc();
2528  refine_uniformly(i_mesh,doc_info);
2529  }
2530 
2531  /// \short p-refine p-refineable sub-meshes, each as many times as
2532  /// specified in the vector and rebuild problem
2533  void p_refine_uniformly(const Vector<unsigned>& nrefine_for_mesh)
2534  {
2535  DocInfo doc_info;
2536  doc_info.disable_doc();
2537  bool prune=false;
2538  p_refine_uniformly_aux(nrefine_for_mesh,doc_info,prune);
2539  }
2540 
2541  /// \short p-refine p-refineable sub-meshes, each as many times as
2542  /// specified in the vector and rebuild problem; doc refinement process
2543  void p_refine_uniformly(const Vector<unsigned>& nrefine_for_mesh,
2544  DocInfo& doc_info)
2545  {
2546  bool prune=false;
2547  p_refine_uniformly_aux(nrefine_for_mesh,doc_info,prune);
2548  }
2549 
2550  /// \short p-refine p-refineable sub-meshes, each as many times as
2551  /// specified in the vector and rebuild problem. Prune after
2552  /// refinements
2553  void p_refine_uniformly_and_prune(const Vector<unsigned>& nrefine_for_mesh)
2554  {
2555  //Not tested:
2556  throw OomphLibWarning("This functionality has not yet been tested.",
2557  "Problem::p_refine_uniformly_and_prune()",
2558  OOMPH_EXCEPTION_LOCATION);
2559  DocInfo doc_info;
2560  doc_info.disable_doc();
2561  bool prune=true;
2562  p_refine_uniformly_aux(nrefine_for_mesh,doc_info,prune);
2563  }
2564 
2565  /// \short p-refine p-refineable sub-meshes, each as many times as
2566  /// specified in the vector and rebuild problem; doc refinement process
2567  void p_refine_uniformly_and_prune(const Vector<unsigned>& nrefine_for_mesh,
2568  DocInfo& doc_info)
2569  {
2570  //Not tested:
2571  throw OomphLibWarning("This functionality has not yet been tested.",
2572  "Problem::p_refine_uniformly_and_prune()",
2573  OOMPH_EXCEPTION_LOCATION);
2574  bool prune=true;
2575  p_refine_uniformly_aux(nrefine_for_mesh,doc_info,prune);
2576  }
2577 
2578  /// \short p-refine (all) p-refineable (sub)mesh(es) uniformly and
2579  /// rebuild problem; doc refinement process.
2580  void p_refine_uniformly(DocInfo& doc_info)
2581  {
2582  // Number of (sub)meshes
2583  unsigned nmesh=std::max(unsigned(1),nsub_mesh());
2584 
2585  // Refine each mesh once
2586  Vector<unsigned> nrefine_for_mesh(nmesh,1);
2587  p_refine_uniformly(nrefine_for_mesh);
2588  }
2589 
2590 
2591  /// \short p-refine (all) p-refineable (sub)mesh(es) uniformly and
2592  /// rebuild problem; doc refinement process.
2594  {
2595  //Not tested:
2596  throw OomphLibWarning("This functionality has not yet been tested.",
2597  "Problem::p_refine_uniformly_and_prune()",
2598  OOMPH_EXCEPTION_LOCATION);
2599  // Number of (sub)meshes
2600  unsigned nmesh=std::max(unsigned(1),nsub_mesh());
2601 
2602  // Refine each mesh once
2603  Vector<unsigned> nrefine_for_mesh(nmesh,1);
2604  p_refine_uniformly_and_prune(nrefine_for_mesh);
2605  }
2606 
2607 
2608  /// \short p-refine (all) p-refineable (sub)mesh(es) uniformly and
2609  /// rebuild problem
2611  {
2612  DocInfo doc_info;
2613  doc_info.disable_doc();
2614  p_refine_uniformly(doc_info);
2615  }
2616 
2617  /// Do uniform p-refinement for submesh i_mesh with documentation
2618  void p_refine_uniformly(const unsigned& i_mesh, DocInfo& doc_info);
2619 
2620  /// Do uniform p-refinement for submesh i_mesh without documentation
2621  void p_refine_uniformly(const unsigned& i_mesh)
2622  {
2623  DocInfo doc_info;
2624  doc_info.disable_doc();
2625  p_refine_uniformly(i_mesh,doc_info);
2626  }
2627 
2628  /// \short Refine (one and only!) mesh by splitting the elements identified
2629  /// by their numbers relative to the problems' only mesh, then rebuild
2630  /// the problem.
2632  elements_to_be_refined);
2633 
2634 
2635  /// \short Refine (one and only!) mesh by splitting the elements identified
2636  /// by their pointers, then rebuild the problem.
2638  elements_to_be_refined_pt);
2639 
2640  /// \short Refine specified submesh by splitting the elements identified
2641  /// by their numbers relative to the submesh, then rebuild the problem.
2642  void refine_selected_elements(const unsigned& i_mesh,
2643  const Vector<unsigned>&
2644  elements_to_be_refined);
2645 
2646  /// \short Refine specified submesh by splitting the elements identified
2647  /// by their pointers, then rebuild the problem.
2648  void refine_selected_elements(const unsigned& i_mesh,
2650  elements_to_be_refined_pt);
2651 
2652  /// \short Refine all submeshes by splitting the elements identified
2653  /// by their numbers relative to each submesh in a Vector of Vectors,
2654  /// then rebuild the problem.
2656  elements_to_be_refined);
2657 
2658  /// \short Refine all submeshes by splitting the elements identified
2659  /// by their pointers within each submesh in a Vector of Vectors,
2660  /// then rebuild the problem.
2662  elements_to_be_refined_pt);
2663 
2664  /// \short p-refine (one and only!) mesh by refining the elements identified
2665  /// by their numbers relative to the problems' only mesh, then rebuild
2666  /// the problem.
2668  elements_to_be_refined);
2669 
2670 
2671  /// \short p-refine (one and only!) mesh by refining the elements identified
2672  /// by their pointers, then rebuild the problem.
2674  elements_to_be_refined_pt);
2675 
2676  /// \short p-refine specified submesh by refining the elements identified
2677  /// by their numbers relative to the submesh, then rebuild the problem.
2678  void p_refine_selected_elements(const unsigned& i_mesh,
2679  const Vector<unsigned>&
2680  elements_to_be_refined);
2681 
2682  /// \short p-refine specified submesh by refining the elements identified
2683  /// by their pointers, then rebuild the problem.
2684  void p_refine_selected_elements(const unsigned& i_mesh,
2686  elements_to_be_refined_pt);
2687 
2688  /// \short p-refine all submeshes by refining the elements identified
2689  /// by their numbers relative to each submesh in a Vector of Vectors,
2690  /// then rebuild the problem.
2692  elements_to_be_refined);
2693 
2694  /// \short p-refine all submeshes by refining the elements identified
2695  /// by their pointers within each submesh in a Vector of Vectors,
2696  /// then rebuild the problem.
2698  elements_to_be_refined_pt);
2699 
2700  /// \short Refine (all) refineable (sub)mesh(es) uniformly and
2701  /// rebuild problem. Return 0 for success,
2702  /// 1 for failure (if unrefinement has reached the coarsest permitted
2703  /// level)
2704  unsigned unrefine_uniformly();
2705 
2706  /// Do uniform refinement for submesh i_mesh without documentation.
2707  /// Return 0 for success, 1 for failure (if unrefinement has reached
2708  /// the coarsest permitted level)
2709  unsigned unrefine_uniformly(const unsigned& i_mesh);
2710 
2711  /// \short p-unrefine (all) p-refineable (sub)mesh(es) uniformly and
2712  /// rebuild problem.
2713  void p_unrefine_uniformly(DocInfo& doc_info);
2714 
2715  /// Do uniform p-unrefinement for submesh i_mesh without documentation.
2716  void p_unrefine_uniformly(const unsigned& i_mesh, DocInfo& doc_info);
2717 
2718  /// \short Adapt problem:
2719  /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2720  /// based on their own error estimates and the target errors specified
2721  /// in the mesh(es). Following mesh adaptation,
2722  /// update global mesh, and re-assign equation numbers.
2723  /// Return # of refined/unrefined elements. On return from this
2724  /// function, Problem can immediately be solved again.
2725  void adapt(unsigned& n_refined, unsigned& n_unrefined);
2726 
2727  /// \short Adapt problem:
2728  /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2729  /// based on their own error estimates and the target errors specified
2730  /// in the mesh(es). Following mesh adaptation,
2731  /// update global mesh, and re-assign equation numbers.
2732  /// On return from this function, Problem can immediately be solved again.
2733  /// [Argument-free wrapper]
2734  void adapt()
2735  {
2736  unsigned n_refined, n_unrefined;
2737  adapt(n_refined,n_unrefined);
2738  }
2739 
2740  /// \short p-adapt problem:
2741  /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2742  /// based on their own error estimates and the target errors specified
2743  /// in the mesh(es). Following mesh adaptation,
2744  /// update global mesh, and re-assign equation numbers.
2745  /// Return # of refined/unrefined elements. On return from this
2746  /// function, Problem can immediately be solved again.
2747  void p_adapt(unsigned& n_refined, unsigned& n_unrefined);
2748 
2749  /// \short p-adapt problem:
2750  /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2751  /// based on their own error estimates and the target errors specified
2752  /// in the mesh(es). Following mesh adaptation,
2753  /// update global mesh, and re-assign equation numbers.
2754  /// On return from this function, Problem can immediately be solved again.
2755  /// [Argument-free wrapper]
2756  void p_adapt()
2757  {
2758  unsigned n_refined, n_unrefined;
2759  p_adapt(n_refined,n_unrefined);
2760  }
2761 
2762 
2763  /// \short Adapt problem:
2764  /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2765  /// based on the error estimates in elemental_error
2766  /// and the target errors specified
2767  /// in the mesh(es). Following mesh adaptation,
2768  /// update global mesh, and re-assign equation numbers.
2769  /// Return # of refined/unrefined elements. On return from this
2770  /// function, Problem can immediately be solved again.
2771  void adapt_based_on_error_estimates(unsigned &n_refined,
2772  unsigned &n_unrefined,
2773  Vector<Vector<double> > &elemental_error);
2774 
2775  /// \short Adapt problem:
2776  /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2777  /// based on the error estimates in elemental_error
2778  /// and the target errors specified
2779  /// in the mesh(es). Following mesh adaptation,
2780  /// update global mesh, and re-assign equation numbers.
2781  /// Return # of refined/unrefined elements. On return from this
2782  /// function, Problem can immediately be solved again.
2783  /// [Wrapper without n_refined and n_unrefined arguments]
2785  Vector<Vector<double> > &elemental_error)
2786  {
2787  unsigned n_refined, n_unrefined;
2788  adapt_based_on_error_estimates(n_refined,n_unrefined,elemental_error);
2789  }
2790 
2791 
2792  /// \short Return the error estimates computed by (all) refineable
2793  /// (sub)mesh(es) in the elemental_error structure, which consists of
2794  /// a vector of vectors of elemental errors, one vector for each (sub)mesh.
2795  void get_all_error_estimates(Vector<Vector<double> > &elemental_error);
2796 
2797  /// \short Get max and min error for all elements in submeshes
2798  void doc_errors(DocInfo& doc_info);
2799 
2800  /// \short Get max and min error for all elements in submeshes
2801  void doc_errors()
2802  {
2803  DocInfo tmp_doc_info;
2804  tmp_doc_info.disable_doc();
2805  doc_errors(tmp_doc_info);
2806  }
2807 
2808  /// \short Enable the output of information when in the newton solver (Default)
2810 
2811  /// \short Disable the output of information when in the newton solver
2813 
2814 };
2815 
2816 
2817 //=======================================================================
2818 /// A class to handle errors in the Newton solver
2819 //=======================================================================
2821 {
2822  public:
2823 
2824  /// Error in the linear solver
2826 
2827  /// Max. # of iterations performed when the Newton solver died
2828  unsigned iterations;
2829 
2830  /// Max. residual when Newton solver died
2831  double maxres;
2832 
2833  /// Default constructor, does nothing
2835  linear_solver_error(false), iterations(0), maxres(0.0) {}
2836 
2837  ///Constructor that passes a failure of the linear solver
2838  NewtonSolverError(const bool& Passed_linear_failure) :
2839  linear_solver_error(Passed_linear_failure), iterations(0), maxres(0.0) {}
2840 
2841  /// Constructor that passes number of iterations and residuals
2842  NewtonSolverError(unsigned Passed_iterations, double Passed_maxres) :
2843  linear_solver_error(false), iterations(Passed_iterations),
2844  maxres(Passed_maxres) {}
2845 
2846 /* /// Broken copy constructor */
2847 /* NewtonSolverError(const NewtonSolverError& dummy) */
2848 /* { */
2849 /* BrokenCopy::broken_copy("NewtonSolverError"); */
2850 /* } */
2851 
2852 /* /// Broken assignment operator */
2853 /* void operator=(const NewtonSolverError&) */
2854 /* { */
2855 /* BrokenCopy::broken_assign("NewtonSolverError"); */
2856 /* } */
2857 
2858 };
2859 
2860 
2861 }
2862 
2863 
2864 #endif
void refine_uniformly(const unsigned &i_mesh)
Do uniform refinement for submesh i_mesh without documentation.
Definition: problem.h:2524
virtual void build_mesh()
Function to build the Problem's base mesh; this must be supplied by the user if they wish to use the ...
Definition: problem.h:1325
virtual void actions_after_explicit_timestep()
Actions that should be performed after each explicit time step.
Definition: problem.h:1064
A class to handle errors in the Newton solver.
Definition: problem.h:2820
void steady_newton_solve(unsigned const &max_adapt=0)
Solve a steady problem using adaptive Newton's method, but in the context of an overall unstady probl...
Definition: problem.cc:9235
Vector< double * > Halo_dof_pt
Storage for the halo degrees of freedom (only required) when accessing via the global equation number...
Definition: problem.h:578
EigenSolver * Default_eigen_solver_pt
Pointer to the default eigensolver.
Definition: problem.h:194
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1234
void remove_duplicate_data(Mesh *const &mesh_pt, bool &actually_removed_some_data)
Private helper function to remove repeated data in external haloed elements in specified mesh...
Definition: problem.cc:2547
Vector< double > Max_res
Maximum residuals at start and after each newton iteration.
Definition: problem.h:603
virtual void actions_after_parameter_increase(double *const &parameter_pt)
Empty virtual function; provides hook to perform actions after the increase in the arclength paramete...
Definition: problem.h:1144
bool is_dparameter_calculated_analytically(double *const &parameter_pt)
Function to determine whether the parameter derivatives are calculated analytically.
Definition: problem.h:280
void newton_solve()
Use Newton method to solve the problem.
Definition: problem.cc:8741
unsigned Sparse_assembly_method
Flag to determine which sparse assembly method to use By default we use assembly by vectors of pairs...
Definition: problem.h:638
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
double Ds_current
Storage for the current step value.
Definition: problem.h:769
virtual void get_inverse_mass_matrix_times_residuals(DoubleVector &Mres)
Return the residual vector multiplied by the inverse mass matrix Virtual so that it can be overloaded...
Definition: problem.cc:3535
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero...
void disable_globally_convergent_newton_method()
disable globally convergent Newton method
Definition: problem.h:1886
Vector< TimeStepper * > Time_stepper_pt
The Vector of time steppers (there could be many different ones in multiphysics problems) ...
Definition: problem.h:204
void delete_all_external_storage()
Wrapper function to delete external storage for each submesh of the problem.
Definition: problem.cc:16221
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 p_refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
p-refine (one and only!) mesh by refining the elements identified by their numbers relative to the pr...
Definition: problem.cc:14989
double Relaxation_factor
Relaxation fator for Newton method (only a fractional Newton correction is applied if this is less th...
Definition: problem.h:589
bool Must_recompute_load_balance_for_assembly
Boolean indicating that the division of elements over processors during the assembly process must be ...
Definition: problem.h:525
void remove_null_pointers_from_external_halo_node_storage()
Consolidate external halo node storage by removing nulled out pointers in external halo and haloed sc...
Definition: problem.cc:3142
void refine_uniformly(const Vector< unsigned > &nrefine_for_mesh)
Refine refineable sub-meshes, each as many times as specified in the vector and rebuild problem...
Definition: problem.h:2448
virtual void actions_before_adapt()
Actions that are to be performed before a mesh adaptation. These might include removing any additiona...
Definition: problem.h:1004
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
virtual ~Problem()
Virtual destructor to clean up memory.
Definition: problem.cc:173
void refine_uniformly_aux(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info, const bool &prune)
Helper function to do compund refinement of (all) refineable (sub)mesh(es) uniformly as many times as...
Definition: problem.cc:15291
void adapt_based_on_error_estimates(unsigned &n_refined, unsigned &n_unrefined, Vector< Vector< double > > &elemental_error)
Adapt problem: Perform mesh adaptation for (all) refineable (sub)mesh(es), based on the error estimat...
Definition: problem.cc:14318
virtual void sparse_assemble_row_or_column_compressed(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Protected helper function that is used to assemble the Jacobian matrix in the case when the storage i...
Definition: problem.cc:4460
LinearSolver * Linear_solver_pt
Pointer to the linear solver for the problem.
Definition: problem.h:176
void doc_errors()
Get max and min error for all elements in submeshes.
Definition: problem.h:2801
void refine_uniformly_and_prune(DocInfo &doc_info)
Refine (all) refineable (sub)mesh(es) uniformly and rebuild problem; doc refinement process...
Definition: problem.h:2500
virtual void sparse_assemble_row_or_column_compressed_with_maps(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
Definition: problem.cc:4567
void refine_uniformly_and_prune(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
Refine refineable sub-meshes, each as many times as specified in the vector and rebuild problem; doc ...
Definition: problem.h:2478
void set_default_first_and_last_element_for_assembly()
Set default first and last elements for parallel assembly of non-distributed problem.
Definition: problem.cc:1597
double arc_length_step_solve_helper(double *const &parameter_pt, const double &ds, const unsigned &max_adapt)
Private helper function that actually contains the guts of the arc-length stepping, parameter_pt is a pointer to the parameter that is traded for the arc-length constraint, ds is the desired arc length and max_adapt is the maximum number of spatial adaptations. The pointer to the parameter may be changed if this is called from the Data-based interface.
Definition: problem.cc:10439
double Newton_solver_tolerance
The Tolerance below which the Newton Method is deemed to have converged.
Definition: problem.h:593
virtual Problem * make_copy()
Make and return a pointer to the copy of the problem. A virtual function that must be filled in by th...
Definition: problem.cc:11857
double dof(const unsigned &i) const
i-th dof in the problem (const version)
Definition: problem.h:1702
virtual double global_temporal_error_norm()
Function to calculate a global error norm, used in adaptive timestepping to control the change in tim...
Definition: problem.h:1202
void disable_store_local_dof_pt_in_elements()
Insist that local dof pointers are NOT set up in each element when equation numbering takes place (th...
Definition: problem.h:1593
double doubly_adaptive_unsteady_newton_solve(const double &dt, const double &epsilon, const unsigned &max_adapt, const bool &first, const bool &shift=true)
Unsteady "doubly" adaptive Newton solve: Does temporal adaptation first, i.e. we try to do a timest...
Definition: problem.h:2347
static ContinuationStorageScheme Continuation_time_stepper
Storage for the single static continuation timestorage object.
Definition: problem.h:750
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 sparse_assemble_row_or_column_compressed_with_lists(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
Definition: problem.cc:4912
const OomphCommunicator * communicator_pt() const
access function to the oomph-lib communicator, const version
Definition: problem.h:1227
friend class BlockFoldLinearSolver
Definition: problem.h:160
virtual void get_eigenproblem_matrices(CRDoubleMatrix &mass_matrix, CRDoubleMatrix &main_matrix, const double &shift=0.0)
Get the matrices required by a eigensolver. If the shift parameter is non-zero the second matrix will...
Definition: problem.cc:8337
bool First_jacobian_sign_change
Boolean to indicate whether a sign change has occured in the Jacobian.
Definition: problem.h:785
double & max_residuals()
Access function to max residuals in Newton iterations before giving up.
Definition: problem.h:1526
void synchronise_dofs(const bool &do_halos, const bool &do_external_halos)
Synchronise the degrees of freedom by overwriting the haloed values with their non-halo counterparts ...
Definition: problem.cc:16328
unsigned iterations
Max. # of iterations performed when the Newton solver died.
Definition: problem.h:2828
void assign_eigenvector_to_dofs(DoubleVector &eigenvector)
Assign the eigenvector passed to the function to the dofs in the problem so that it can be output by ...
Definition: problem.cc:8671
Information for documentation of results: Directory and file number to enable output in the form RESL...
double doubly_adaptive_unsteady_newton_solve(const double &dt, const double &epsilon, const unsigned &max_adapt, const unsigned &suppress_resolve_after_spatial_adapt_flag, const bool &first, const bool &shift=true)
Unsteady "doubly" adaptive Newton solve: Does temporal adaptation first, i.e. we try to do a timest...
Definition: problem.h:2380
EigenSolver *& eigen_solver_pt()
Return a pointer to the eigen solver object.
Definition: problem.h:1440
std::map< double *, bool > Calculate_dparameter_analytic
Map used to determine whether the derivatives with respect to a parameter should be finite difference...
Definition: problem.h:243
void enable_mass_matrix_reuse()
Enable recycling of the mass matrix in explicit timestepping schemes. Useful for timestepping on fixe...
Definition: problem.cc:11648
virtual void actions_after_read_unstructured_meshes()
Actions that are to be performed before reading in restart data for problems involving unstructured b...
Definition: problem.h:1092
Mesh *const & mesh_pt(const unsigned &imesh) const
Return a pointer to the i-th submesh (const version)
Definition: problem.h:1280
bool Bypass_increase_in_dof_check_during_pruning
Boolean to bypass check of increase in dofs during pruning.
Definition: problem.h:954
cstr elem_len * i
Definition: cfortran.h:607
void p_refine_uniformly_and_prune(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
p-refine p-refineable sub-meshes, each as many times as specified in the vector and rebuild problem; ...
Definition: problem.h:2567
unsigned long set_timestepper_for_all_data(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data=false)
Set all problem data to have the same timestepper (timestepper_pt) Return the new number of dofs in t...
Definition: problem.cc:11412
NewtonSolverError(unsigned Passed_iterations, double Passed_maxres)
Constructor that passes number of iterations and residuals.
Definition: problem.h:2842
double adaptive_unsteady_newton_solve(const double &dt_desired, const double &epsilon)
Attempt to advance timestep by dt_desired. If the solution fails the timestep will be halved until co...
Definition: problem.cc:10941
Class to keep track of discrete/continous time. It is essential to have a single Time object when usi...
Definition: timesteppers.h:67
void copy(Problem *orig_problem_pt)
Copy Data values, nodal positions etc from specified problem. Note: This is not a copy constructor...
Definition: problem.cc:11706
double DTSF_max_increase
Maximum possible increase of dt between time-steps in adaptive schemes.
Definition: problem.h:706
void set_consistent_pinned_values_for_continuation()
Private helper function that is used to set the appropriate pinned values for continuation.
Definition: problem.cc:10389
Vector< Mesh * > Sub_mesh_pt
Vector of pointers to submeshes.
Definition: problem.h:173
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned...
virtual void actions_after_adapt()
Actions that are to be performed after a mesh adaptation.
Definition: problem.h:1007
bool Arc_length_step_taken
Boolean to indicate whether an arc-length step has been taken.
Definition: problem.h:788
void flush_global_data()
Flush the Problem's global data – resizes container to zero. Data objects are not deleted! ...
Definition: problem.h:1567
void unset_analytic_dparameter(double *const &parameter_pt)
Function to turn off analytic calculation of the parameter derivatives in continuation and bifurcatio...
Definition: problem.h:266
void load_balance(DocInfo &doc_info, const bool &report_stats)
Balance the load of a (possibly non-uniformly refined) problem that has already been distributed...
Definition: problem.h:1383
double Max_permitted_error_for_halo_check
Threshold for error throwing in Problem::check_halo_schemes()
Definition: problem.h:2053
bool Discontinuous_element_formulation
Is the problem a discontinuous one, i.e. can the elemental contributions be treated independently...
Definition: problem.h:693
The Problem class.
Definition: problem.h:152
void describe_dofs(std::ostream &out=*(oomph_info.stream_pt())) const
Function to describe the dofs in terms of the global equation number, i.e. what type of value (nodal ...
Definition: problem.cc:2334
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
Vector< Data * > Global_data_pt
Vector of global data: "Nobody" (i.e. none of the elements etc.) is "in charge" of this Data so it wo...
Definition: problem.h:423
double *& dof_pt(const unsigned &i)
Pointer to i-th dof in the problem.
Definition: problem.h:1708
void set_explicit_time_stepper_pt(ExplicitTimeStepper *const &explicit_time_stepper_pt)
Set the explicit timestepper for the problem. The function will automatically create or resize the Ti...
Definition: problem.cc:1569
Vector< double > Dof_current
Storage for the present values of the variables.
Definition: problem.h:766
void refine_distributed_base_mesh(Vector< Vector< Vector< unsigned > > > &to_be_refined_on_each_root, const unsigned &max_level_overall)
Load balance helper routine: refine each new base (sub)mesh based upon the elements to be refined wit...
Definition: problem.cc:19747
bool Jacobian_reuse_is_enabled
Is re-use of Jacobian in Newton iteration enabled? Default: false.
Definition: problem.h:615
bool Shut_up_in_newton_solve
Boolean to indicate if all output is suppressed in Problem::newton_solve(). Defaults to false...
Definition: problem.h:2069
char t
Definition: cfortran.h:572
virtual void actions_after_implicit_timestep_and_error_estimation()
Actions that should be performed after each implicit time step. This is needed if your actions_after_...
Definition: problem.h:1058
void(* SpatialErrorEstimatorWithDocFctPt)(Mesh *&mesh_pt, Vector< double > &elemental_error, DocInfo &doc_info)
Function pointer for spatial error estimator with doc.
Definition: problem.h:1237
bool Calculate_hessian_products_analytic
Map used to determine whether the hessian products should be computed using finite differences...
Definition: problem.h:248
Vector< GeneralisedElement * > Base_mesh_element_pt
Vector to store the correspondence between a root element and its element number within the global me...
Definition: problem.h:544
Vector< double > elemental_assembly_time()
Return vector of most-recent elemental assembly times (used for load balancing). Zero sized if no Jac...
Definition: problem.h:852
bool Use_globally_convergent_newton_method
Use the globally convergent newton method.
Definition: problem.h:217
virtual void read(std::ifstream &restart_file, bool &unsteady_restart)
Read refinement pattern of all refineable meshes and refine them accordingly, then read all Data and ...
Definition: problem.cc:12100
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
double Max_residuals
Maximum desired residual: if the maximum residual exceeds this value, the program will exit...
Definition: problem.h:607
void(* SpatialErrorEstimatorFctPt)(Mesh *&mesh_pt, Vector< double > &elemental_error)
Function pointer for spatial error estimator.
Definition: problem.h:1233
void get_all_halo_data(std::map< unsigned, double * > &map_of_halo_data)
Get pointers to all possible halo data indexed by global equation number in a map.
Definition: problem.cc:16243
void enable_discontinuous_formulation()
Indicate that the problem involves discontinuous elements This allows for a more efficiently assembly...
Definition: problem.h:1615
virtual void sparse_assemble_row_or_column_compressed_with_vectors_of_pairs(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
Definition: problem.cc:5327
bool Empty_actions_after_read_unstructured_meshes_has_been_called
Boolean to indicate that empty actions_after_read_unstructured_meshes() function has been called...
Definition: problem.h:225
double * bifurcation_parameter_pt() const
Return pointer to the parameter that is used in the bifurcation detection. If we are not tracking a b...
Definition: problem.cc:9991
double Parameter_current
Storage for the present value of the global parameter.
Definition: problem.h:744
bool Empty_actions_before_read_unstructured_meshes_has_been_called
Boolean to indicate that empty actions_before_read_unstructured_meshes() function has been called...
Definition: problem.h:221
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1264
void p_refine_uniformly(DocInfo &doc_info)
p-refine (all) p-refineable (sub)mesh(es) uniformly and rebuild problem; doc refinement process...
Definition: problem.h:2580
void p_refine_uniformly_aux(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info, const bool &prune)
Helper function to do compund p-refinement of (all) p-refineable (sub)mesh(es) uniformly as many time...
Definition: problem.cc:15442
OomphInfo oomph_info
Distributed_problem_matrix_distribution & distributed_problem_matrix_distribution()
accesss function to the distributed matrix distribution method 1 - Automatic - the Problem distributi...
Definition: problem.h:834
void disable_discontinuous_formulation()
Disable the use of a discontinuous-element formulation. Note that the methods will all still work eve...
Definition: problem.h:1622
void calculate_predictions()
Calculate predictions.
Definition: problem.cc:11496
void copy_haloed_eqn_numbers_helper(const bool &do_halos, const bool &do_external_halos)
A private helper function to copy the haloed equation numbers into the halo equation numbers...
Definition: problem.cc:16791
void set_dofs(const DoubleVector &dofs)
Set the values of the dofs.
Definition: problem.cc:3361
unsigned & max_newton_iterations()
Access function to max Newton iterations before giving up.
Definition: problem.h:1515
void clear_elemental_assembly_time()
Clear storage of most-recent elemental assembly times (used for load balancing). Next load balancing ...
Definition: problem.h:858
void synchronise_all_dofs()
Perform all required synchronisation in solvers.
Definition: problem.cc:16305
Distributed_problem_matrix_distribution Dist_problem_matrix_distribution
The distributed matrix distribution method 1 - Automatic - the Problem distribution is employed...
Definition: problem.h:942
unsigned unrefine_uniformly()
Refine (all) refineable (sub)mesh(es) uniformly and rebuild problem. Return 0 for success...
Definition: problem.cc:15693
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1221
virtual void actions_after_change_in_global_parameter(double *const &parameter_pt)
Actions that are to be performed when the global parameter addressed by parameter_pt has been changed...
Definition: problem.h:1116
double Minimum_dt
Minimum desired dt: if dt falls below this value, exit.
Definition: problem.h:699
virtual void shift_time_values()
Shift all values along to prepare for next timestep.
Definition: problem.cc:11473
void refine_uniformly(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
Refine refineable sub-meshes, each as many times as specified in the vector and rebuild problem; doc ...
Definition: problem.h:2458
void recompute_load_balanced_assembly()
Helper function to re-assign the first and last elements to be assembled by each processor during par...
Definition: problem.cc:1676
LinearSolver *& mass_matrix_solver_for_explicit_timestepper_pt()
Definition: problem.h:1431
bool & time_adaptive_newton_crash_on_solve_fail()
Access function for Time_adaptive_newton_crash_on_solve_fail.
Definition: problem.h:1529
void problem_is_nonlinear(const bool &prob_lin)
Access function to Problem_is_nonlinear.
Definition: problem.h:1518
double & minimum_dt()
Access function to min timestep in adaptive timestepping.
Definition: problem.h:1509
double * dof_pt(const unsigned &i) const
Pointer to i-th dof in the problem (const version)
Definition: problem.h:1711
virtual void debug_hook_fct(const unsigned &i)
Hook for debugging. Can be overloaded in driver code; argument allows identification of where we're c...
Definition: problem.h:254
bool Use_default_partition_in_load_balance
Flag to use "default partition" during load balance. Should only be set to true when run in validatio...
Definition: problem.h:510
virtual void symmetrise_eigenfunction_for_adaptive_pitchfork_tracking()
Virtual function that is used to symmetrise the problem so that the current solution exactly satisfie...
Definition: problem.cc:9969
ExplicitTimeStepper * Explicit_time_stepper_pt
Pointer to a single explicit timestepper.
Definition: problem.h:207
void enable_doc_imbalance_in_parallel_assembly()
Enable doc of load imbalance in parallel assembly of distributed problem.
Definition: problem.h:841
virtual void read(std::ifstream &restart_file)
Read refinement pattern of all refineable meshes and refine them accordingly, then read all Data and ...
Definition: problem.h:1944
void p_refine_uniformly(const Vector< unsigned > &nrefine_for_mesh)
p-refine p-refineable sub-meshes, each as many times as specified in the vector and rebuild problem ...
Definition: problem.h:2533
NewtonSolverError(const bool &Passed_linear_failure)
Constructor that passes a failure of the linear solver.
Definition: problem.h:2838
void disable_info_in_newton_solve()
Disable the output of information when in the newton solver.
Definition: problem.h:2812
void p_unrefine_uniformly(DocInfo &doc_info)
p-unrefine (all) p-refineable (sub)mesh(es) uniformly and rebuild problem.
Definition: problem.cc:15827
bool Mass_matrix_reuse_is_enabled
Is re-use of the mass matrix in explicit timestepping enabled Default:false.
Definition: problem.h:684
bool & use_predictor_values_as_initial_guess()
Definition: problem.h:1871
unsigned ntime_stepper() const
Return the number of time steppers.
Definition: problem.h:1578
virtual void sparse_assemble_row_or_column_compressed_with_two_vectors(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
Definition: problem.cc:5686
void set_analytic_dparameter(double *const &parameter_pt)
Function to turn on analytic calculation of the parameter derivatives in continuation and bifurcation...
Definition: problem.h:261
void add_time_stepper_pt(TimeStepper *const &time_stepper_pt)
Add a timestepper to the problem. The function will automatically create or resize the Time object so...
Definition: problem.cc:1528
void get_dofs(DoubleVector &dofs) const
Return the vector of dofs, i.e. a vector containing the current values of all unknowns.
Definition: problem.cc:2454
Mesh * Mesh_pt
The mesh pointer.
Definition: problem.h:170
void refine_uniformly(DocInfo &doc_info)
Refine (all) refineable (sub)mesh(es) uniformly and rebuild problem; doc refinement process...
Definition: problem.h:2487
bool Bisect_to_find_bifurcation
Boolean to control wheter bisection is used to located bifurcation.
Definition: problem.h:782
double Continuation_direction
The direction of the change in parameter that will ensure that a branch is followed in one direction ...
Definition: problem.h:738
bool linear_solver_error
Error in the linear solver.
Definition: problem.h:2825
AssemblyHandler *& assembly_handler_pt()
Return a pointer to the assembly handler object.
Definition: problem.h:1502
Time * Time_pt
Pointer to global time for the problem.
Definition: problem.h:200
void calculate_continuation_derivatives_helper(const DoubleVector &z)
A function that performs the guts of the continuation derivative calculation in arc length continuati...
Definition: problem.cc:9834
unsigned setup_element_count_per_dof()
Function that populates the Element_counter_per_dof vector with the number of elements that contribut...
Definition: problem.cc:218
void disable_jacobian_reuse()
Disable recycling of Jacobian in Newton iteration.
Definition: problem.h:1859
virtual void get_dvaluesdt(DoubleVector &f)
Get the time derivative of all values (using get_inverse_mass_matrix_times_residuals(..) with all time steppers set to steady) e.g. for use in explicit time steps. The approach used is slighty hacky, beware if you have a residual which is non-linear or implicit in the derivative or if you have overloaded get_jacobian(...).
Definition: problem.cc:3640
double Numerical_zero_for_sparse_assembly
A tolerance used to determine whether the entry in a sparse matrix is zero. If it is then storage nee...
Definition: problem.h:665
virtual void get_jacobian(DoubleVector &residuals, SumOfMatrices &jacobian)
Dummy virtual function that must be overloaded by the problem to specify which matrices should be sum...
Definition: problem.h:1755
EigenSolver * Eigen_solver_pt
Pointer to the eigen solver for the problem.
Definition: problem.h:185
double Desired_proportion_of_arc_length
Proportion of the arc-length to taken by the parameter.
Definition: problem.h:727
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
unsigned Dof_current_offset
Definition: problem.h:760
void refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
Refine (one and only!) mesh by splitting the elements identified by their numbers relative to the pro...
Definition: problem.cc:14710
LinearAlgebraDistribution *const & dof_distribution_pt() const
Return the pointer to the dof distribution (read-only)
Definition: problem.h:1570
void assign_initial_values_impulsive()
Initialise data and nodal positions to simulate impulsive start from initial configuration/solution.
Definition: problem.cc:11346
double & dof(const unsigned &i)
i-th dof in the problem
Definition: problem.h:1696
void disable_doc_imbalance_in_parallel_assembly()
Disable doc of load imbalance in parallel assembly of distributed problem.
Definition: problem.h:846
double & dof_derivative(const unsigned &i)
Access function to the derivative of the i-th (local) dof with respect to the arc length...
Definition: problem.h:1150
double doubly_adaptive_unsteady_newton_solve_helper(const double &dt, const double &epsilon, const unsigned &max_adapt, const unsigned &suppress_resolve_after_spatial_adapt, const bool &first, const bool &shift=true)
Private helper function that actually performs the unsteady "doubly" adaptive Newton solve...
Definition: problem.cc:11230
void prune_halo_elements_and_nodes(DocInfo &doc_info, const bool &report_stats)
(Irreversibly) prune halo(ed) elements and nodes, usually after another round of refinement, to get rid of excessively wide halo layers. Note that the current mesh will be now regarded as the base mesh and no unrefinement relative to it will be possible once this function has been called.
Definition: problem.cc:1032
double & newton_solver_tolerance()
Access function to tolererance of the Newton solver, i.e. the maximum value of the residuals that wil...
Definition: problem.h:1534
OomphCommunicator * Communicator_pt
The communicator for this problem.
Definition: problem.h:1216
void store_current_dof_values()
Store the current values of the degrees of freedom.
Definition: problem.cc:8567
unsigned Sparse_assemble_with_arrays_initial_allocation
the number of elements to initially allocate for a matrix row within the sparse_assembly_with_two_arr...
Definition: problem.h:652
std::ostream *& stream_pt()
Access function for the stream pointer.
void enable_problem_distributed()
Enable problem distributed.
Definition: problem.h:959
void set_pinned_values_to_zero()
Set all pinned values to zero. Used to set boundary conditions to be homogeneous in the copy of the p...
Definition: problem.cc:4294
void calculate_continuation_derivatives_fd(double *const &parameter_pt)
A function to calculate the derivatives with respect to the arc-length required for continuation by f...
Definition: problem.cc:9748
Time * time_pt() const
Return a pointer to the global time object (const version).
Definition: problem.h:1449
void explicit_timestep(const double &dt, const bool &shift_values=true)
Take an explicit timestep of size dt and optionally shift any stored values of the time history...
Definition: problem.cc:10814
double Theta_squared
Value of the scaling parameter required so that the parameter occupies the desired proportion of the ...
Definition: problem.h:731
void p_refine_uniformly_and_prune(const Vector< unsigned > &nrefine_for_mesh)
p-refine p-refineable sub-meshes, each as many times as specified in the vector and rebuild problem...
Definition: problem.h:2553
void prune_halo_elements_and_nodes(const bool &report_stats=false)
(Irreversibly) prune halo(ed) elements and nodes, usually after another round of refinement, to get rid of excessively wide halo layers. Note that the current mesh will be now regarded as the base mesh and no unrefinement relative to it will be possible once this function has been called.
Definition: problem.h:2045
double Parameter_derivative
Storage for the derivative of the global parameter wrt arc-length.
Definition: problem.h:741
bool Doc_time_in_distribute
Protected boolean flag to provide comprehensive timimings during problem distribution. Initialised to false.
Definition: problem.h:634
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
bool Time_adaptive_newton_crash_on_solve_fail
Bool to specify what to do if a Newton solve fails within a time adaptive solve. Default (false) is t...
Definition: problem.h:612
void reset_assembly_handler_to_default()
Reset the system to the standard non-augemented state.
Definition: problem.cc:10191
void get_all_error_estimates(Vector< Vector< double > > &elemental_error)
Return the error estimates computed by (all) refineable (sub)mesh(es) in the elemental_error structur...
Definition: problem.cc:14421
void dump(const std::string &dump_file_name) const
Dump refinement pattern of all refineable meshes and all generic Problem data to file for restart...
Definition: problem.h:1956
double Maximum_dt
Maximum desired dt.
Definition: problem.h:702
unsigned Sparse_assemble_with_arrays_allocation_increment
the number of elements to add to a matrix row when the initial allocation is exceeded within the spar...
Definition: problem.h:657
virtual void actions_before_explicit_timestep()
Actions that should be performed before each explicit time step.
Definition: problem.h:1061
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1574
bool Mass_matrix_has_been_computed
Has the mass matrix been computed (and can therefore be reused) Default: false.
Definition: problem.h:688
virtual void actions_after_distribute()
Actions to be performed after a (mesh) distribution.
Definition: problem.h:1102
void p_refine_uniformly(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
p-refine p-refineable sub-meshes, each as many times as specified in the vector and rebuild problem; ...
Definition: problem.h:2543
unsigned Nnewton_iter_taken
Actual number of Newton iterations taken during the most recent iteration.
Definition: problem.h:600
void get_flat_packed_refinement_pattern_for_load_balancing(const Vector< unsigned > &old_domain_for_base_element, const Vector< unsigned > &new_domain_for_base_element, const unsigned &max_refinement_level_overall, std::map< unsigned, Vector< unsigned > > &flat_packed_refinement_info_for_root)
Get flat-packed refinement pattern for each root element in current mesh (labeled by unique number of...
Definition: problem.cc:19595
std::map< GeneralisedElement *, unsigned > Base_mesh_element_number_plus_one
Map which stores the correspondence between a root element and its element number (plus one) within t...
Definition: problem.h:534
unsigned nrow() const
access function to the number of global rows.
void reset_arc_length_parameters()
Reset the "internal" arc-length continuation parameters, so as to allow continuation in another param...
Definition: problem.h:2286
DoubleVectorWithHaloEntries Element_count_per_dof
Counter that records how many elements contribute to each dof. Used to determine the (discrete) arc-l...
Definition: problem.h:557
void solve_eigenproblem(const unsigned &n_eval, Vector< std::complex< double > > &eigenvalue, const bool &steady=true)
Solve an eigenproblem as assembled by EigenElements, but only return the eigenvalues, not the eigenvectors. The boolean flag (default true) is used to specify whether the weighted mass-matrix terms from the timestepping scheme should be included in the jacobian.
Definition: problem.h:1814
double & maximum_dt()
Access function to max timestep in adaptive timestepping.
Definition: problem.h:1512
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1446
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void activate_pitchfork_tracking(double *const &parameter_pt, const DoubleVector &symmetry_vector, const bool &block_solve=true)
Turn on pitchfork tracking using the augmented system specified in the PitchForkHandler class...
Definition: problem.cc:10102
void activate_fold_tracking(double *const &parameter_pt, const bool &block_solve=true)
Turn on fold tracking using the augmented system specified in the FoldHandler class. After a call to this function subsequent calls of the standard solution methods will converge to a fold (limit) point at a particular value of the variable addressed by parameter_pt. The system may not converge if the initial guess is sufficiently poor or, alternatively, if finite differencing is used to calculate the jacobian matrix in the elements. If the boolean flag block_solver is true (the default) then a block factorisation is used to solve the augmented system which is both faster and uses less memory.
Definition: problem.cc:10011
NewtonSolverError()
Default constructor, does nothing.
Definition: problem.h:2834
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
TimeStepper *& time_stepper_pt()
Access function for the pointer to the first (presumably only) timestepper.
Definition: problem.h:1460
bool Pause_at_end_of_sparse_assembly
Protected boolean flag to halt program execution during sparse assemble process to assess peak memory...
Definition: problem.h:630
void setup_dof_halo_scheme()
Function that is used to setup the halo scheme.
Definition: problem.cc:284
void enable_info_in_newton_solve()
Enable the output of information when in the newton solver (Default)
Definition: problem.h:2809
Vector< unsigned > First_el_for_assembly
First element to be assembled by given processor for non-distributed problem (only kept up to date wh...
Definition: problem.h:515
void get_fd_jacobian(DoubleVector &residuals, DenseMatrix< double > &jacobian)
Return the fully-assembled Jacobian and residuals, generated by finite differences.
Definition: problem.cc:7804
void set_first_and_last_element_for_assembly(Vector< unsigned > &first_el_for_assembly, Vector< unsigned > &last_el_for_assembly)
Manually set first and last elements for parallel assembly of non-distributed problem.
Definition: problem.h:970
Distributed_problem_matrix_distribution
enum for distribution of distributed jacobians. 1 - Automatic - the Problem distribution is employed...
Definition: problem.h:819
LinearSolver * Default_linear_solver_pt
Pointer to the default linear solver.
Definition: problem.h:191
void p_adapt()
p-adapt problem: Perform mesh adaptation for (all) refineable (sub)mesh(es), based on their own error...
Definition: problem.h:2756
void p_refine_uniformly(const unsigned &i_mesh)
Do uniform p-refinement for submesh i_mesh without documentation.
Definition: problem.h:2621
Vector< double > Elemental_assembly_time
Storage for assembly times (used for load balancing)
Definition: problem.h:570
void p_refine_uniformly()
p-refine (all) p-refineable (sub)mesh(es) uniformly and rebuild problem
Definition: problem.h:2610
int Sign_of_jacobian
Storage for the sign of the global Jacobian.
Definition: problem.h:734
unsigned nglobal_data() const
Return the number of global data values.
Definition: problem.h:1581
A class that is used to define the functions used to assemble the elemental contributions to the resi...
virtual void sparse_assemble_row_or_column_compressed_with_two_arrays(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
Definition: problem.cc:6052
double & dof_current(const unsigned &i)
Access function to the current value of the i-th (local) dof at the start of a continuation step...
Definition: problem.h:1160
static bool Suppress_warning_about_actions_before_read_unstructured_meshes
Flag to allow suppression of warning messages re reading in unstructured meshes during restart...
Definition: problem.h:313
bool Use_continuation_timestepper
Boolean to control original or new storage of dof stuff.
Definition: problem.h:747
virtual void actions_before_newton_solve()
Any actions that are to be performed before a complete Newton solve (e.g. adjust boundary conditions)...
Definition: problem.h:1015
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
void disable_mass_matrix_reuse()
Turn off recyling of the mass matrix in explicit timestepping schemes.
Definition: problem.cc:11673
Vector< double > * Saved_dof_pt
Pointer to vector for backup of dofs.
Definition: problem.h:210
GeneralisedTimestepper used to store the arclength derivatives and pervious solutions required in con...
void get_my_eqns(AssemblyHandler *const &assembly_handler_pt, const unsigned &el_lo, const unsigned &el_hi, Vector< unsigned > &my_eqns)
Helper method that returns the (unique) global equations to which the elements in the range el_lo to ...
Definition: problem.cc:6504
void get_data_to_be_sent_during_load_balancing(const Vector< unsigned > &element_domain_on_this_proc, Vector< int > &send_n, Vector< double > &send_data, Vector< int > &send_displacement, Vector< unsigned > &old_domain_for_base_element, Vector< unsigned > &new_domain_for_base_element, unsigned &max_refinement_level_overall)
Load balance helper routine: Get data to be sent to other processors during load balancing and other ...
Definition: problem.cc:19048
bool does_pointer_correspond_to_problem_data(double *const &parameter_pt)
Return a boolean flag to indicate whether the pointer parameter_pt refers to values stored in a Data ...
Definition: problem.cc:9777
double Timestep_reduction_factor_after_nonconvergence
What it says: If temporally adaptive Newton solver fails to to converge, reduce timestep by this fact...
Definition: problem.h:2081
void check_halo_schemes()
Check the halo/haloed node/element schemes.
Definition: problem.h:1995
void initialise_dt(const double &dt)
Set all timesteps to the same value, dt, and assign weights for all timesteppers in the problem...
Definition: problem.cc:13001
double FD_step_used_in_get_hessian_vector_products
Definition: problem.h:679
unsigned local_index(const unsigned &global_eqn)
Return the local index associated with the global equation.
unsigned add_sub_mesh(Mesh *const &mesh_pt)
Add a submesh to the problem and return its number, i, by which it can be accessed via mesh_pt(i)...
Definition: problem.h:1293
void unset_default_partition_in_load_balance()
Do not use the default partition in the load balance.
Definition: problem.h:1411
A class for compressed column matrices that store doubles.
Definition: matrices.h:2573
void calculate_continuation_derivatives_fd_helper(double *const &parameter_pt)
A function that performs the guts of the continuation derivative calculation in arc-length continuati...
Definition: problem.cc:9917
virtual void get_jacobian(DoubleVector &residuals, DenseDoubleMatrix &jacobian)
Return the fully-assembled Jacobian and residuals for the problem Interface for the case when the Jac...
Definition: problem.cc:3851
LinearSolver * mass_matrix_solver_for_explicit_timestepper_pt() const
Definition: problem.h:1436
void disable_problem_distributed()
Disable problem distributed.
Definition: problem.h:962
void add_eigenvector_to_dofs(const double &epsilon, const DoubleVector &eigenvector)
Add the eigenvector passed to the function scaled by the constat epsilon to the dofs in the problem s...
Definition: problem.cc:8704
void disable_doc()
Disable documentation.
double maxres
Max. residual when Newton solver died.
Definition: problem.h:2831
void flush_sub_meshes()
Flush the problem's collection of sub-meshes. Must be followed by call to rebuild_global_mesh().
Definition: problem.h:1302
Problem()
Constructor: Allocate space for one time stepper and set all pointers to NULL and set defaults for al...
Definition: problem.cc:75
double arc_length_step_solve(double *const &parameter_pt, const double &ds, const unsigned &max_adapt=0)
Solve a steady problem using arc-length continuation, when the parameter that becomes a variable corr...
Definition: problem.cc:10210
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: problem.h:1424
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
bool jacobian_reuse_is_enabled()
Is recycling of Jacobian in Newton iteration enabled?
Definition: problem.h:1866
void send_data_to_be_sent_during_load_balancing(Vector< int > &send_n, Vector< double > &send_data, Vector< int > &send_displacement)
Load balance helper routine: Send data to other processors during load balancing. ...
Definition: problem.cc:18726
Mesh *& mesh_pt(const unsigned &imesh)
Return a pointer to the i-th submesh. If there are no submeshes, the 0-th submesh is the global mesh ...
Definition: problem.h:1271
bool Bifurcation_detection
Boolean to control bifurcation detection via determinant of Jacobian.
Definition: problem.h:779
bool Use_finite_differences_for_continuation_derivatives
Definition: problem.h:792
Assembly_method
Enumerated flags to determine which sparse assembly method is used.
Definition: problem.h:641
void refine_uniformly_and_prune(const Vector< unsigned > &nrefine_for_mesh)
Refine refineable sub-meshes, each as many times as specified in the vector and rebuild problem...
Definition: problem.h:2468
Problem(const Problem &dummy)
Broken copy constructor.
Definition: problem.h:1248
void operator=(const Problem &)
Broken assignment operator.
Definition: problem.h:1254
bool Default_set_initial_condition_called
Has default set_initial_condition function been called? Default: false.
Definition: problem.h:214
bool Jacobian_has_been_computed
Has a Jacobian been computed (and can therefore be re-used if required)? Default: false...
Definition: problem.h:619
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
Vector< double > Dof_derivative
Storage for the derivative of the problem variables wrt arc-length.
Definition: problem.h:763
unsigned Max_newton_iterations
Maximum number of Newton iterations.
Definition: problem.h:596
void load_balance(const bool &report_stats)
Balance the load of a (possibly non-uniformly refined) problem that has already been distributed...
Definition: problem.h:1364
unsigned nsub_mesh() const
Return number of submeshes.
Definition: problem.h:1289
Data *& global_data_pt(const unsigned &i)
Return a pointer to the the i-th global data object.
Definition: problem.h:1557
void unsteady_newton_solve(const double &dt)
Advance time by dt and solve by Newton's method. This version always shifts time values.
Definition: problem.cc:10843
void bifurcation_adapt_helper(unsigned &n_refined, unsigned &n_unrefined, const unsigned &bifurcation_type, const bool &actually_adapt=true)
A function that is used to adapt a bifurcation-tracking problem, which requires separate interpolatio...
Definition: problem.cc:13119
virtual void partition_global_mesh(Mesh *&global_mesh_pt, DocInfo &doc_info, Vector< unsigned > &element_domain, const bool &report_stats=false)
Definition: problem.cc:863
Vector< unsigned > distribute(const Vector< unsigned > &element_partition, DocInfo &doc_info, const bool &report_stats=false)
Distribute the problem and doc, using the specified partition; returns a vector which details the par...
Definition: problem.cc:390
bool Use_predictor_values_as_initial_guess
Use values from the time stepper predictor as an initial guess.
Definition: problem.h:232
void p_refine_uniformly_and_prune(DocInfo &doc_info)
p-refine (all) p-refineable (sub)mesh(es) uniformly and rebuild problem; doc refinement process...
Definition: problem.h:2593
virtual void dump(std::ofstream &dump_file) const
Dump refinement pattern of all refineable meshes and all generic Problem data to file for restart...
Definition: problem.cc:11876
unsigned newton_solve_continuation(double *const &parameter_pt)
Perform a basic arc-length continuation step using Newton's method. Returns number of Newton steps ta...
Definition: problem.cc:9323
void refine_uniformly()
Refine (all) refineable (sub)mesh(es) uniformly and rebuild problem.
Definition: problem.h:2513
void build_global_mesh()
Build the global mesh by combining the all the submeshes. Note: The nodes boundary information refers...
Definition: problem.cc:1472
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
Definition: problem.cc:1965
virtual void actions_before_newton_step()
Any actions that are to be performed before each individual Newton step. Most likely to be used for d...
Definition: problem.h:1036
AssemblyHandler * Default_assembly_handler_pt
Pointer to the default assembly handler.
Definition: problem.h:197
bool Always_take_one_newton_step
Boolean to indicate whether a Newton step should be taken even if the initial residuals are below the...
Definition: problem.h:2076
void setup_base_mesh_info_after_pruning()
Helper function to re-setup the Base_mesh enumeration (used during load balancing) after pruning...
Definition: problem.cc:19860
Vector< double * > Dof_pt
Vector of pointers to dofs.
Definition: problem.h:551
void load_balance()
Balance the load of a (possibly non-uniformly refined) problem that has already been distributed...
Definition: problem.h:1343
DoubleVectorHaloScheme * Halo_scheme_pt
Pointer to the halo scheme for any global vectors that have the Dof_distribution. ...
Definition: problem.h:574
virtual void set_initial_condition()
Set initial condition (incl previous timesteps). We need to establish this interface because I...
Definition: problem.h:1172
void get_first_and_last_element_for_assembly(Vector< unsigned > &first_el_for_assembly, Vector< unsigned > &last_el_for_assembly) const
Get first and last elements for parallel assembly of non-distributed problem.
Definition: problem.h:986
bool problem_has_been_distributed()
Access to Problem_has_been_distributed flag.
Definition: problem.h:2056
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
double DTSF_min_decrease
Minimum allowed decrease of dt between time-steps in adaptive schemes. Lower scaling values will reje...
Definition: problem.h:711
bool mass_matrix_reuse_is_enabled()
Return whether the mass matrix is being reused.
Definition: problem.h:2444
EigenSolver *const & eigen_solver_pt() const
Return a pointer to the eigen solver object (const version)
Definition: problem.h:1443
virtual void actions_before_read_unstructured_meshes()
Actions that are to be performed before reading in restart data for problems involving unstructured b...
Definition: problem.h:1078
unsigned self_test()
Self-test: Check meshes and global data. Return 0 for OK.
Definition: problem.cc:13046
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
void adapt_based_on_error_estimates(Vector< Vector< double > > &elemental_error)
Adapt problem: Perform mesh adaptation for (all) refineable (sub)mesh(es), based on the error estimat...
Definition: problem.h:2784
virtual void actions_before_distribute()
Actions to be performed before a (mesh) distribution.
Definition: problem.h:1099
bool are_hessian_products_calculated_analytically()
Function to determine whether the hessian products are calculated analytically.
Definition: problem.h:301
LinearSolver * Mass_matrix_solver_for_explicit_timestepper_pt
Definition: problem.h:182
bool Problem_is_nonlinear
Boolean flag indicating if we're dealing with a linear or nonlinear Problem – if set to false the New...
Definition: problem.h:625
void get_bifurcation_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction calculated as part of a bifurcation tracking process. If we are not tracking...
Definition: problem.cc:9999
long synchronise_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Classify any non-classified nodes into halo/haloed and synchronise equation numbers. Return the total number of degrees of freedom in the overall problem.
Definition: problem.cc:16554
const TimeStepper * time_stepper_pt() const
Access function for the pointer to the first (presumably only) timestepper.
Definition: problem.h:1473
bool Store_local_dof_pt_in_elements
Boolean to indicate whether local dof pointers should be stored in the elements.
Definition: problem.h:229
virtual void actions_after_newton_step()
Any actions that are to be performed after each individual Newton step. Most likely to be used for di...
Definition: problem.h:1041
AssemblyHandler * Assembly_handler_pt
Definition: problem.h:188
void activate_bifurcation_tracking(double *const &parameter_pt, const DoubleVector &eigenvector, const bool &block_solve=true)
Activate generic bifurcation tracking for a single (real) eigenvalue where the initial guess for the ...
Definition: problem.cc:10037
unsigned Desired_newton_iterations_ds
The desired number of Newton Steps to reach convergence at each step along the arc.
Definition: problem.h:773
void calculate_continuation_derivatives(double *const &parameter_pt)
A function to calculate the derivatives wrt the arc-length. This version of the function actually doe...
Definition: problem.cc:9644
LinearSolver *const & linear_solver_pt() const
Return a pointer to the linear solver object (const version)
Definition: problem.h:1427
void activate_hopf_tracking(double *const &parameter_pt, const bool &block_solve=true)
Turn on Hopf bifurcation tracking using the augmented system specified in the HopfHandler class...
Definition: problem.cc:10135
ExplicitTimeStepper *& explicit_time_stepper_pt()
Return a pointer to the explicit timestepper.
Definition: problem.h:1489
A Base class for explicit timesteppers.
void deactivate_bifurcation_tracking()
Deactivate all bifuraction tracking, by reseting the assembly handler to the default.
Definition: problem.h:2228
void unset_analytic_hessian_products()
Function to turn off analytic calculation of the parameter derivatives in continuation and bifurcatio...
Definition: problem.h:296
virtual void actions_after_newton_solve()
Any actions that are to be performed after a complete Newton solve, e.g. post processing. CAREFUL: This step should (and if the FD-based LinearSolver FD_LU is used, must) only update values that are pinned!
Definition: problem.h:1021
Vector< Problem * > Copy_of_problem_pt
Vector of pointers to copies of the problem used in adaptive bifurcation tracking problems (ALH: TEMP...
Definition: problem.h:238
unsigned Parallel_sparse_assemble_previous_allocation
The amount of data allocated during the previous parallel sparse assemble. Used to optimise the next ...
Definition: problem.h:946
void solve_eigenproblem(const unsigned &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector, const bool &steady=true)
Get derivative of an element in the problem wrt a global parameter, used in continuation problems...
Definition: problem.cc:8282
Mesh *const & mesh_pt() const
Return a pointer to the global mesh (const version)
Definition: problem.h:1267
void parallel_sparse_assemble(const LinearAlgebraDistribution *const &dist_pt, Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residuals)
Helper method to assemble CRDoubleMatrices from distributed on multiple processors.
Definition: problem.cc:6551
bool distributed() const
Definition: problem.h:798
void enable_store_local_dof_pt_in_elements()
Insist that local dof pointers are set up in each element when equation numbering takes place...
Definition: problem.h:1588
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
Vector< unsigned > Last_el_plus_one_for_assembly
Last element (plus one) to be assembled by given processor for non-distributed problem (only kept up ...
Definition: problem.h:520
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
void add_global_data(Data *const &global_data_pt)
Add Data to the Problem's global data – the Problem will perform equation numbering etc...
Definition: problem.h:1561
bool Problem_has_been_distributed
Has the problem been distributed amongst multiple processors?
Definition: problem.h:951
LinearAlgebraDistribution *& distribution_pt()
Return the pointer to the distirbution used to setup the halo information.
void globally_convergent_line_search(const Vector< double > &x_old, const double &half_residual_squared_old, DoubleVector &gradient, DoubleVector &newton_dir, double &half_residual_squared, const double &stpmax)
Line search helper for globally convergent Newton method.
Definition: problem.cc:9095
virtual void actions_after_implicit_timestep()
Actions that should be performed after each implicit time step. This is needed when one wants to solv...
Definition: problem.h:1053
void enable_jacobian_reuse()
Enable recycling of Jacobian in Newton iteration (if the linear solver allows it). Useful for linear problems with constant Jacobians or nonlinear problems where you're willing to risk the trade-off between faster solve times and degraded Newton convergence rate.
Definition: problem.h:1852
void send_refinement_info_helper(Vector< unsigned > &old_domain_for_base_element, Vector< unsigned > &new_domain_for_base_element, const unsigned &max_refinement_level_overall, std::map< unsigned, Vector< unsigned > > &refinement_info_for_root_local, Vector< Vector< Vector< unsigned > > > &refinement_info_for_root_elements)
Send refinement information between processors.
Definition: problem.cc:18134
bool Doc_imbalance_in_parallel_assembly
Boolean to switch on assessment of load imbalance in parallel assembly of distributed problem...
Definition: problem.h:506
void set_default_partition_in_load_balance()
Set the use of the default partition in the load balance.
Definition: problem.h:1407
A general mesh class.
Definition: mesh.h:74
void enable_globally_convergent_newton_method()
enable globally convergent Newton method
Definition: problem.h:1880
void restore_dof_values()
Restore the stored values of the degrees of freedom.
Definition: problem.cc:8609
void adapt()
Adapt problem: Perform mesh adaptation for (all) refineable (sub)mesh(es), based on their own error e...
Definition: problem.h:2734
double Minimum_dt_but_still_proceed
If Minimum_dt_but_still_proceed positive, then dt will not be reduced below this value during adaptiv...
Definition: problem.h:718
void set_analytic_hessian_products()
Function to turn on analytic calculation of the parameter derivatives in continuation and bifurcation...
Definition: problem.h:291
void rebuild_global_mesh()
If one of the submeshes has changed (e.g. by mesh adaptation) we need to update the global mesh...
Definition: problem.cc:1516
double & time()
Return the current value of continuous time.
Definition: problem.cc:11377
TimeStepper *& time_stepper_pt(const unsigned &i)
Return a pointer to the i-th timestepper.
Definition: problem.h:1485
virtual void actions_before_implicit_timestep()
Actions that should be performed before each implicit time step. This is needed when one wants to sol...
Definition: problem.h:1047
bool Scale_arc_length
Boolean to control whether arc-length should be scaled.
Definition: problem.h:724
double Minimum_ds
Minimum desired value of arc-length.
Definition: problem.h:776
void add_to_dofs(const double &lambda, const DoubleVector &increment_dofs)
Add lambda x incremenet_dofs[l] to the l-th dof.
Definition: problem.cc:3519
unsigned Dof_derivative_offset
Definition: problem.h:755
AssemblyHandler *const & assembly_handler_pt() const
Return a pointer to the assembly handler object (const version)
Definition: problem.h:1505
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:57
void bifurcation_adapt_doc_errors(const unsigned &bifurcation_type)
A function that is used to document the errors used in the adaptive solution of bifurcation problems...
Definition: problem.cc:13410