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