problem.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1302 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-10-23 12:52:19 +0100 (Mon, 23 Oct 2017) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
31 #ifdef OOMPH_HAS_MPI
32 #include "mpi.h"
33 #endif
34 
35 #include<list>
36 #include<algorithm>
37 #include<string>
38 
39 #include "oomph_utilities.h"
40 #include "problem.h"
41 #include "timesteppers.h"
42 #include "explicit_timesteppers.h"
44 #include "refineable_mesh.h"
45 #include "triangle_mesh.h"
46 #include "linear_solver.h"
47 #include "eigen_solver.h"
48 #include "assembly_handler.h"
49 #include "dg_elements.h"
50 #include "partitioning.h"
51 #include "spines.h"
52 
53 //Include to fill in additional_setup_shared_node_scheme() function
55 
56 
57 namespace oomph
58 {
59 
60 
61 //////////////////////////////////////////////////////////////////
62 //Non-inline functions for the problem class
63 //////////////////////////////////////////////////////////////////
64 
65 //=================================================================
66 ///The continuation timestepper object
67 //=================================================================
68 ContinuationStorageScheme Problem::Continuation_time_stepper;
69 
70 //================================================================
71 /// Constructor: Allocate space for one time stepper
72 /// and set all pointers to NULL and set defaults for all
73 /// parameters.
74 //===============================================================
76  Mesh_pt(0), Time_pt(0), Explicit_time_stepper_pt(0), Saved_dof_pt(0),
77  Default_set_initial_condition_called(false),
78  Use_globally_convergent_newton_method(false),
79  Empty_actions_before_read_unstructured_meshes_has_been_called(false),
80  Empty_actions_after_read_unstructured_meshes_has_been_called(false),
81  Store_local_dof_pt_in_elements(false),
82  Calculate_hessian_products_analytic(false),
83 #ifdef OOMPH_HAS_MPI
84  Doc_imbalance_in_parallel_assembly(false),
85  Use_default_partition_in_load_balance(false),
86  Must_recompute_load_balance_for_assembly(true),
87  Halo_scheme_pt(0),
88 #endif
89  Relaxation_factor(1.0),
90  Newton_solver_tolerance(1.0e-8),
92  Nnewton_iter_taken(0),
93  Max_residuals(10.0),
94  Time_adaptive_newton_crash_on_solve_fail(false),
95  Jacobian_reuse_is_enabled(false), Jacobian_has_been_computed(false),
96  Problem_is_nonlinear(true),
97  Pause_at_end_of_sparse_assembly(false),
98  Doc_time_in_distribute(false),
99  Sparse_assembly_method(Perform_assembly_using_vectors_of_pairs),
100  Sparse_assemble_with_arrays_initial_allocation(400),
101  Sparse_assemble_with_arrays_allocation_increment(150),
102  Numerical_zero_for_sparse_assembly(0.0),
103  FD_step_used_in_get_hessian_vector_products(1.0e-8),
104  Mass_matrix_reuse_is_enabled(false), Mass_matrix_has_been_computed(false),
105  Discontinuous_element_formulation(false),
106  Minimum_dt(1.0e-12), Maximum_dt(1.0e12),
107  DTSF_max_increase(4.0),
108  DTSF_min_decrease(0.8),
109  Minimum_dt_but_still_proceed(-1.0),
110  Scale_arc_length(true), Desired_proportion_of_arc_length(0.5),
111  Theta_squared(1.0), Sign_of_jacobian(0), Continuation_direction(1.0),
112  Parameter_derivative(1.0), Parameter_current(0.0),
113  Use_continuation_timestepper(false),
114  Dof_derivative_offset(1),
115  Dof_current_offset(2),
116  Ds_current(0.0), Desired_newton_iterations_ds(5),
117  Minimum_ds(1.0e-10), Bifurcation_detection(false),
118  Bisect_to_find_bifurcation(false),
119  First_jacobian_sign_change(false),Arc_length_step_taken(false),
120  Use_finite_differences_for_continuation_derivatives(false),
121 #ifdef OOMPH_HAS_MPI
122  Dist_problem_matrix_distribution(Uniform_matrix_distribution),
123  Parallel_sparse_assemble_previous_allocation(0),
124  Problem_has_been_distributed(false),
125  Bypass_increase_in_dof_check_during_pruning(false),
126  Max_permitted_error_for_halo_check(1.0e-14),
127 #endif
128  Shut_up_in_newton_solve(false),
129  Always_take_one_newton_step(false),
130  Timestep_reduction_factor_after_nonconvergence(0.5)
131  {
133 
134  /// Setup terminate helper
136 
137  // By default no submeshes:
138  Sub_mesh_pt.resize(0);
139  // No timesteppers
140  Time_stepper_pt.resize(0);
141 
142  //Set the linear solvers, eigensolver and assembly handler
145 
147 
149 
150  // setup the communicator
151 #ifdef OOMPH_HAS_MPI
153  {
155  }
156  else
157  {
159  }
160 #else
162 #endif
163 
164  // just create an empty linear algebra distribution for the
165  // DOFs
166  // this is setup when assign_eqn_numbers(...) is called.
168  }
169 
170 //================================================================
171 /// Destructor to clean up memory
172 //================================================================
174  {
175 
176  // Delete the memory assigned for the global time
177  // (it's created on the fly in Problem::add_time_stepper_pt()
178  // so we are entitled to delete it.
179  if (Time_pt!=0)
180  {
181  delete Time_pt;
182  Time_pt = 0;
183  }
184 
185  // We're not using the default linear solver,
186  // somebody else must have built it, so that person
187  // must be in charge of killing it.
188  // We can safely delete the defaults, however
190 
193  delete Communicator_pt;
194  delete Dof_distribution_pt;
195 
196  // Delete any copies of the problem that have been created for
197  // use in adaptive bifurcation tracking.
198  // ALH: This will eventually go
199  unsigned n_copies = Copy_of_problem_pt.size();
200  for(unsigned c=0;c<n_copies;c++)
201  {
202  delete Copy_of_problem_pt[c];
203  }
204 
205  // if this problem has sub meshes then we must delete the Mesh_pt
206  if (Sub_mesh_pt.size() != 0)
207  {
209  delete Mesh_pt;
210  }
211  }
212 
213  //=================================================================
214  ///Setup the count vector that records how many elements contribute
215  ///to each degree of freedom. Returns the total number of elements
216  ///in the problem
217  //=================================================================
219  {
220  //Now set the element counter to have the current Dof distribution
222  //We need to use the halo scheme (assuming it has been setup)
223 #ifdef OOMPH_HAS_MPI
225 #endif
226 
227  //Loop over the elements and count the entries
228  //and number of (non-halo) elements
229  const unsigned n_element = this->mesh_pt()->nelement();
230  unsigned n_non_halo_element_local=0;
231  for(unsigned e=0;e<n_element;e++)
232  {
233  GeneralisedElement* elem_pt = this->mesh_pt()->element_pt(e);
234 #ifdef OOMPH_HAS_MPI
235  //Ignore halo elements
236  if(!elem_pt->is_halo())
237  {
238 #endif
239  //Increment the number of non halo elements
240  ++n_non_halo_element_local;
241  //Now count the number of times the element contributes to a value
242  //using the current assembly handler
243  unsigned n_var = this->Assembly_handler_pt->ndof(elem_pt);
244  for(unsigned n=0;n<n_var;n++)
245  {
247  this->Assembly_handler_pt->eqn_number(elem_pt,n));
248  }
249 #ifdef OOMPH_HAS_MPI
250  }
251 #endif
252  }
253 
254  //Storage for the total number of elements
255  unsigned Nelement=0;
256 
257  //Add together all the counts if we are in parallel
258 #ifdef OOMPH_HAS_MPI
260 
261  //If distributed, find the total number of elements in the problem
263  {
264  //Need to gather the total number of non halo elements
265  MPI_Allreduce(&n_non_halo_element_local,&Nelement,1,MPI_UNSIGNED,MPI_SUM,
266  this->communicator_pt()->mpi_comm());
267  }
268  //Otherwise the total number is the same on each processor
269  else
270 #endif
271  {
272  Nelement = n_non_halo_element_local;
273  }
274 
275  return Nelement;
276  }
277 
278 
279 #ifdef OOMPH_HAS_MPI
280 
281  //==================================================================
282  /// Setup the halo scheme for the degrees of freedom
283  //==================================================================
285  {
286  //Find the number of elements stored on this processor
287  const unsigned n_element = this->mesh_pt()->nelement();
288 
289  //Work out the all global equations to which this processor
290  //contributes
291  Vector<unsigned> my_eqns;
292  this->get_my_eqns(this->Assembly_handler_pt,0,n_element-1,my_eqns);
293 
294  //Build the halo scheme, based on the equations to which this
295  //processor contributes
297  new DoubleVectorHaloScheme(this->Dof_distribution_pt,my_eqns);
298 
299  //Find pointers to all the halo dofs
300  //There may be more of these than required by my_eqns
301  //(but should not be less)
302  std::map<unsigned,double*> halo_data_pt;
303  this->get_all_halo_data(halo_data_pt);
304 
305  //Now setup the Halo_dofs
306  Halo_scheme_pt->setup_halo_dofs(halo_data_pt,this->Halo_dof_pt);
307  }
308 
309  //==================================================================
310  /// Distribute the problem without doc; report stats if required.
311  /// Returns actual partitioning used, e.g. for restart.
312  //==================================================================
313  Vector<unsigned> Problem::distribute(const bool& report_stats)
314  {
315  // Set dummy doc paramemters
316  DocInfo doc_info;
317  doc_info.disable_doc();
318 
319  // Set the sizes of the input and output vectors
320  unsigned n_element=mesh_pt()->nelement();
321  Vector<unsigned> element_partition(n_element,0);
322 
323  // Distribute and return partitioning
324  return distribute(element_partition,doc_info,report_stats);
325  }
326 
327  //==================================================================
328  /// Distribute the problem according to specified partition.
329  /// If all entries in partitioning vector are zero we use METIS
330  /// to do the partitioning after all.
331  /// Returns actual partitioning used, e.g. for restart.
332  //==================================================================
334  (const Vector<unsigned>& element_partition, const bool& report_stats)
335  {
336 #ifdef PARANOID
337  bool has_non_zero_entry=false;
338  unsigned n=element_partition.size();
339  for (unsigned i=0;i<n;i++)
340  {
341  if (element_partition[i]!=0)
342  {
343  has_non_zero_entry=true;
344  break;
345  }
346  }
347  if (!has_non_zero_entry)
348  {
349  std::ostringstream warn_message;
350  warn_message << "WARNING: All entries in specified partitioning vector \n"
351  << " are zero -- will ignore this and use METIS\n"
352  << " to perform the partitioning\n";
353  OomphLibWarning(warn_message.str(),
354  "Problem::distribute()",
355  OOMPH_EXCEPTION_LOCATION);
356  }
357 #endif
358  // Set dummy doc paramemters
359  DocInfo doc_info;
360  doc_info.disable_doc();
361 
362  // Distribute and return partitioning
363  return distribute(element_partition,doc_info,report_stats);
364  }
365 
366  //==================================================================
367  /// Distribute the problem and doc to specified DocInfo.
368  /// Returns actual partitioning used, e.g. for restart.
369  //==================================================================
371  const bool& report_stats)
372  {
373  // Set the sizes of the input and output vectors
374  unsigned n_element=mesh_pt()->nelement();
375 
376  // Dummy input vector
377  Vector<unsigned> element_partition(n_element,0);
378 
379  // Distribute and return partitioning
380  return distribute(element_partition,doc_info,report_stats);
381  }
382 
383  //==================================================================
384  /// Distribute the problem according to specified partition.
385  /// (If all entries in partitioning vector are zero we use METIS
386  /// to do the partitioning after all) and doc.
387  /// Returns actual partitioning used, e.g. for restart.
388  //==================================================================
390  (const Vector<unsigned>& element_partition,
391  DocInfo& doc_info, const bool& report_stats)
392  {
393  // Storage for number of processors and number of elements in global mesh
394  int n_proc=this->communicator_pt()->nproc();
395  int my_rank=this->communicator_pt()->my_rank();
396  int n_element=mesh_pt()->nelement();
397 
398  // Vector to be returned
399  Vector<unsigned> return_element_domain;
400 
401  // Buffer extreme cases
402  if (n_proc==1) // single-process job - don't do anything
403  {
404  if (report_stats)
405  {
406  std::ostringstream warn_message;
407  warn_message << "WARNING: You've tried to distribute a problem over\n"
408  << "only one processor: this would make METIS crash.\n"
409  << "Ignoring your request for distribution.\n";
410  OomphLibWarning(warn_message.str(),
411  "Problem::distribute()",
412  OOMPH_EXCEPTION_LOCATION);
413  }
414  }
415  else if (n_proc>n_element) // more processors than elements
416  {
417  // Throw an error
418  std::ostringstream error_stream;
419  error_stream << "You have tried to distribute a problem\n"
420  << "but there are less elements than processors.\n"
421  << "Please re-run with more elements!\n"
422  << "Please also ensure that actions_before_distribute().\n"
423  << "and actions_after_distribute() are correctly set up.\n"
424  << std::endl;
425  throw OomphLibError(error_stream.str(),
426  OOMPH_CURRENT_FUNCTION,
427  OOMPH_EXCEPTION_LOCATION);
428  }
429  else
430  {
431  // We only distribute uniformly-refined meshes; buffer the case where
432  // either mesh is not uniformly refined
433  bool a_mesh_is_not_uniformly_refined=false;
434  unsigned n_mesh=nsub_mesh();
435  if (n_mesh==0)
436  {
437  // Check refinement levels
438  if (TreeBasedRefineableMeshBase* mmesh_pt =
439  dynamic_cast<TreeBasedRefineableMeshBase*>(mesh_pt(0)))
440  {
441  unsigned min_ref_level=0;
442  unsigned max_ref_level=0;
443  mmesh_pt->get_refinement_levels(min_ref_level,max_ref_level);
444  // If they are not the same
445  if (max_ref_level!=min_ref_level)
446  {
447  a_mesh_is_not_uniformly_refined=true;
448  }
449  }
450  }
451  else
452  {
453  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
454  {
455  // Check refinement levels for each mesh individually
456  // (one mesh is allowed to be "more uniformly refined" than another)
457  if (TreeBasedRefineableMeshBase* mmesh_pt =
458  dynamic_cast<TreeBasedRefineableMeshBase*>(mesh_pt(i_mesh)))
459  {
460  unsigned min_ref_level=0;
461  unsigned max_ref_level=0;
462  mmesh_pt->get_refinement_levels(min_ref_level,max_ref_level);
463  // If they are not the same
464  if (max_ref_level!=min_ref_level)
465  {
466  a_mesh_is_not_uniformly_refined=true;
467  }
468  }
469  }
470  }
471 
472  // If any mesh is not uniformly refined
473  if (a_mesh_is_not_uniformly_refined)
474  {
475  // Again it may make more sense to throw an error here as the user
476  // will probably not be running a problem that is small enough to
477  // fit the whole of on each processor
478  std::ostringstream error_stream;
479  error_stream << "You have tried to distribute a problem\n"
480  << "but at least one of your meshes is no longer\n"
481  << "uniformly refined. In order to preserve the Tree\n"
482  << "and TreeForest structure, Problem::distribute() can\n"
483  << "only be called while meshes are uniformly refined.\n"
484  << std::endl;
485  throw OomphLibError(error_stream.str(),
486  OOMPH_CURRENT_FUNCTION,
487  OOMPH_EXCEPTION_LOCATION);
488  }
489  else
490  {
491  // Is there any global data? If so, distributing the problem won't work
492  if (nglobal_data() > 0)
493  {
494  std::ostringstream error_stream;
495  error_stream << "You have tried to distribute a problem\n"
496  << "and there is some global data.\n"
497  << "This is not likely to work...\n"
498  << std::endl;
499  throw OomphLibError(error_stream.str(),
500  OOMPH_CURRENT_FUNCTION,
501  OOMPH_EXCEPTION_LOCATION);
502  }
503 
504  double t_start = 0;
505  if (Doc_time_in_distribute)
506  {
507  t_start=TimingHelpers::timer();
508  }
509 
510 
511 #ifdef PARANOID
512  unsigned old_ndof=ndof();
513 #endif
514 
515  // Need to partition the global mesh before distributing
516  Mesh* global_mesh_pt = mesh_pt();
517 
518  // Vector listing the affiliation of each element
519  unsigned nelem=global_mesh_pt->nelement();
520  Vector<unsigned> element_domain(nelem);
521 
522  // Number of elements that I'm in charge of, based on any
523  // incoming partitioning
524  unsigned n_my_elements=0;
525 
526  // Have we used the pre-set partitioning
527  bool used_preset_partitioning=false;
528 
529  // Partition the mesh, unless the partition has already been passed in
530  // If it hasn't then the sum of all the entries of the vector should be 0
531  unsigned sum_element_partition=0;
532  unsigned n_part=element_partition.size();
533  for (unsigned e=0;e<n_part;e++)
534  {
535  // ... another one for me.
536  if (int(element_partition[e])==my_rank) n_my_elements++;
537 
538  sum_element_partition+=element_partition[e];
539  }
540  if (sum_element_partition==0)
541  {
542  oomph_info << "INFO: using METIS to partition elements"
543  << std::endl;
544  partition_global_mesh(global_mesh_pt,doc_info,element_domain);
545  used_preset_partitioning=false;
546  }
547  else
548  {
549  oomph_info << "INFO: using pre-set partition of elements"
550  << std::endl;
551  used_preset_partitioning=true;
552  element_domain=element_partition;
553  }
554 
555  // Set the GLOBAL Mesh as being distributed
556  global_mesh_pt->set_communicator_pt(this->communicator_pt());
557 
558  double t_end = 0.0;
559  if (Doc_time_in_distribute)
560  {
561  t_end=TimingHelpers::timer();
562  oomph_info << "Time for partitioning of global mesh: "
563  << t_end-t_start << std::endl;
564  t_start = TimingHelpers::timer();
565  }
566 
567  // Store how many elements we had in the various sub-meshes
568  // before actions_before_distribute() (which may empty some of
569  // them).
570  Vector<unsigned> n_element_in_old_submesh(n_mesh);
571  if (n_mesh!=0)
572  {
573  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
574  {
575  unsigned nsub_elem=mesh_pt(i_mesh)->nelement();
576  n_element_in_old_submesh[i_mesh]=nsub_elem;
577  }
578  }
579 
580  // Partitioning complete; call actions before distribute
581  actions_before_distribute();
582 
583  if (Doc_time_in_distribute)
584  {
585  t_end = TimingHelpers::timer();
586  oomph_info << "Time for actions before distribute: "
587  << t_end-t_start << std::endl;
588  }
589 
590  // This next bit is cheap -- omit timing
591  // t_start = TimingHelpers::timer();
592 
593  // Number of submeshes (NB: some may have been deleted in
594  // actions_after_distribute())
595  n_mesh=nsub_mesh();
596 
597 
598  // Prepare vector of vectors for submesh element domains
599  Vector<Vector<unsigned> > submesh_element_domain(n_mesh);
600 
601  // The submeshes need to know their own element domains.
602  // Also if any meshes have been emptied we ignore their
603  // partitioning in the vector that we return from here
604  return_element_domain.reserve(element_domain.size());
605  if (n_mesh!=0)
606  {
607  unsigned count=0;
608  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
609  {
610  unsigned nsub_elem=mesh_pt(i_mesh)->nelement();
611  submesh_element_domain[i_mesh].resize(nsub_elem);
612  unsigned nsub_elem_old=n_element_in_old_submesh[i_mesh];
613  for (unsigned e=0; e<nsub_elem_old; e++)
614  {
615  if (nsub_elem_old==nsub_elem)
616  {
617  submesh_element_domain[i_mesh][e]=element_domain[count];
618  return_element_domain.push_back(element_domain[count]);
619  }
620  //return_element_domain.push_back(element_domain[count]);
621  count++;
622  }
623  }
624  }
625  else
626  {
627  return_element_domain=element_domain;
628  }
629 
630  if (Doc_time_in_distribute)
631  {
632  t_start = TimingHelpers::timer();
633  }
634 
635  // Setup the map between "root" element and number in global mesh
636  // (currently used in the load_balance() routines)
637 
638  // This map is only established for structured meshes, then we
639  // need to check here the type of mesh
640  if (n_mesh==0)
641  {
642  // Check if the only one mesh is an structured mesh
643  bool structured_mesh = true;
644  TriangleMeshBase* tri_mesh_pt =
645  dynamic_cast<TriangleMeshBase*>(mesh_pt(0));
646  if (tri_mesh_pt != 0)
647  {
648  structured_mesh = false;
649  } // if (tri_mesh_pt != 0)
650  if (structured_mesh)
651  {
652  const unsigned n_ele = global_mesh_pt->nelement();
653  Base_mesh_element_pt.resize(n_ele);
654  Base_mesh_element_number_plus_one.clear();
655  for (unsigned e=0;e<n_ele;e++)
656  {
657  GeneralisedElement* el_pt=global_mesh_pt->element_pt(e);
658  Base_mesh_element_number_plus_one[el_pt]=e+1;
659  Base_mesh_element_pt[e]=el_pt;
660  } // for (e<n_ele)
661  } // A TreeBaseMesh mesh
662  } // if (n_mesh==0)
663  else
664  {
665  // If we have submeshes then we only add those elements that
666  // belong to structured meshes, but first compute the number
667  // of total elements in the structured meshes
668  unsigned nglobal_element = 0;
669  // Store which submeshes are structured
670  std::vector<bool> is_structured_mesh(n_mesh);
671  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
672  {
673  TriangleMeshBase* tri_mesh_pt =
674  dynamic_cast<TriangleMeshBase*>(mesh_pt(i_mesh));
675  if (tri_mesh_pt != 0)
676  {
677  // Set the flags to indicate this is not an structured
678  // mesh
679  is_structured_mesh[i_mesh] = false;
680  } // if (tri_mesh_pt != 0)
681  else
682  {
683  // Set the flags to indicate this is an structured
684  // mesh
685  is_structured_mesh[i_mesh] = true;
686  } // else if (tri_mesh_pt != 0)
687  // Check if mesh is an structured mesh
688  if (is_structured_mesh[i_mesh])
689  {
690  nglobal_element+=mesh_pt(i_mesh)->nelement();
691  } // A TreeBaseMesh mesh
692  } // for (i_mesh<n_mesh)
693 
694  // Once computed the number of elements, then resize the
695  // structure
696  Base_mesh_element_pt.resize(nglobal_element);
697  Base_mesh_element_number_plus_one.clear();
698  unsigned counter = 0;
699  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
700  {
701  // Check if mesh is an structured mesh
702  if (is_structured_mesh[i_mesh])
703  {
704  const unsigned n_ele = mesh_pt(i_mesh)->nelement();
705  for (unsigned e=0;e<n_ele;e++)
706  {
707  GeneralisedElement* el_pt=mesh_pt(i_mesh)->element_pt(e);
708  Base_mesh_element_number_plus_one[el_pt]=counter+1;
709  Base_mesh_element_pt[counter]=el_pt;
710  // Inrease the global element number
711  counter++;
712  } // for (e<n_ele)
713  } // An structured mesh
714  } // for (i_mesh<n_mesh)
715 
716 #ifdef PARANOID
717  if (counter != nglobal_element)
718  {
719  std::ostringstream error_stream;
720  error_stream
721  << "The number of global elements ("<<nglobal_element
722  <<") is not the sameas the number of\nadded elements ("
723  << counter <<") to the Base_mesh_element_pt data "
724  << "structure!!!\n\n";
725  throw OomphLibError(error_stream.str(),
726  "Problem::distribute()",
727  OOMPH_EXCEPTION_LOCATION);
728  } // if (counter != nglobal_element)
729 #endif // #ifdef PARANOID
730 
731  } // else if (n_mesh==0)
732 
733  // Wipe everything if a pre-determined partitioning
734  // didn't specify ANY elements for this processor
735  // (typically happens during restarts with larger number
736  // of processors -- in this case we really want an empty
737  // processor rather than one with any "kept" halo elements)
738  bool overrule_keep_as_halo_element_status=false;
739  if ((n_my_elements==0)&&(used_preset_partitioning))
740  {
741  oomph_info << "INFO: We're over-ruling the \"keep as halo element\"\n"
742  << " status because the preset partitioning\n"
743  << " didn't place ANY elements on this processor,\n"
744  << " probably because of a restart on a larger \n"
745  << " number of processors\n";
746  overrule_keep_as_halo_element_status=true;
747  }
748 
749 
750  // Distribute the (sub)meshes (i.e. sort out their halo lookup schemes)
751  Vector<GeneralisedElement*> deleted_element_pt;
752  if (n_mesh==0)
753  {
754  global_mesh_pt->distribute(this->communicator_pt(),
755  element_domain,
756  deleted_element_pt,
757  doc_info,report_stats,
758  overrule_keep_as_halo_element_status);
759  }
760  else // There are submeshes, "distribute" each one separately
761  {
762  for (unsigned i_mesh=0; i_mesh<n_mesh; i_mesh++)
763  {
764  if (report_stats)
765  {
766  oomph_info << "Distributing submesh " << i_mesh << std::endl
767  << "--------------------" << std::endl;
768  }
769  // Set the doc_info number to reflect the submesh
770  doc_info.number()=i_mesh;
771  mesh_pt(i_mesh)->distribute(this->communicator_pt(),
772  submesh_element_domain[i_mesh],
773  deleted_element_pt,
774  doc_info,report_stats,
775  overrule_keep_as_halo_element_status);
776  }
777  // Rebuild the global mesh
778  rebuild_global_mesh();
779  }
780 
781  // Null out information associated with deleted elements
782  unsigned n_del=deleted_element_pt.size();
783  for (unsigned e=0;e<n_del;e++)
784  {
785  GeneralisedElement* el_pt=deleted_element_pt[e];
786  unsigned old_el_number=Base_mesh_element_number_plus_one[el_pt]-1;
787  Base_mesh_element_number_plus_one[el_pt]=0;
788  Base_mesh_element_pt[old_el_number]=0;
789  }
790 
791  if (Doc_time_in_distribute)
792  {
793  t_end = TimingHelpers::timer();
794  oomph_info << "Time for mesh-level distribution: "
795  << t_end-t_start << std::endl;
796  t_start = TimingHelpers::timer();
797  }
798 
799  // Now the problem has been distributed
800  Problem_has_been_distributed=true;
801 
802  // Call actions after distribute
803  actions_after_distribute();
804 
805  if (Doc_time_in_distribute)
806  {
807  t_end = TimingHelpers::timer();
808  oomph_info << "Time for actions after distribute: "
809  << t_end-t_start << std::endl;
810  t_start = TimingHelpers::timer();
811  }
812 
813  // Re-assign the equation numbers (incl synchronisation if reqd)
814  unsigned n_dof=assign_eqn_numbers();
815  oomph_info << "Number of equations: " << n_dof
816  << std::endl;
817 
818  if (Doc_time_in_distribute)
819  {
820  t_end = TimingHelpers::timer();
821  oomph_info << "Time for re-assigning eqn numbers (in distribute): "
822  << t_end-t_start << std::endl;
823  }
824 
825 
826 #ifdef PARANOID
827  if (n_dof!=old_ndof)
828  {
829  std::ostringstream error_stream;
830  error_stream
831  << "Number of dofs in distribute() has changed "
832  << "from " << old_ndof << " to " << n_dof << "\n"
833  <<"Check that you've implemented any necessary actions_before/after\n"
834  << "distribute functions, e.g. to pin redundant pressure dofs"
835  << " etc.\n";
836  throw OomphLibError(error_stream.str(),
837  OOMPH_CURRENT_FUNCTION,
838  OOMPH_EXCEPTION_LOCATION);
839  }
840 #endif
841 
842  } // end if to check for uniformly refined mesh(es)
843 
844  } // end if to check number of processors vs. number of elements etc.
845 
846 
847  // Force re-analysis of time spent on assembly each
848  // elemental Jacobian
849  Must_recompute_load_balance_for_assembly=true;
850  Elemental_assembly_time.clear();
851 
852  // Return the partition vector used in the distribution
853  return return_element_domain;
854 
855  }
856 
857  //==================================================================
858  /// Partition the global mesh, return vector specifying the processor
859  /// number for each element. Virtual so that it can be overloaded by
860  /// any user; the default is to use METIS to perform the partitioning
861  /// (with a bit of cleaning up afterwards to sort out "special cases").
862  //==================================================================
863  void Problem::partition_global_mesh(Mesh* &global_mesh_pt, DocInfo& doc_info,
864  Vector<unsigned>& element_domain,
865  const bool& report_stats)
866  {
867  // Storage for number of processors and current processor
868  int n_proc=this->communicator_pt()->nproc();
869  int rank=this->communicator_pt()->my_rank();
870 
871  std::ostringstream filename;
872  std::ofstream some_file;
873 
874  // Doc the original mesh on proc 0
875  //--------------------------------
876  if (doc_info.is_doc_enabled())
877  {
878  if (rank==0)
879  {
880  filename << doc_info.directory() << "/complete_mesh"
881  << doc_info.number() << ".dat";
882  global_mesh_pt->output(filename.str().c_str(),5);
883  }
884  }
885 
886  // Partition the mesh
887  //-------------------
888  // METIS Objective (0: minimise edge cut; 1: minimise total comm volume)
889  unsigned objective=0;
890 
891  // Do the partitioning
892  unsigned nelem=0;
893  if (this->communicator_pt()->my_rank()==0)
894  {
895  METIS::partition_mesh(this,n_proc,objective,element_domain);
896  nelem=element_domain.size();
897  }
898  MPI_Bcast(&nelem,1,MPI_UNSIGNED,0,this->communicator_pt()->mpi_comm());
899  element_domain.resize(nelem);
900  MPI_Bcast(&element_domain[0],nelem,MPI_UNSIGNED,0,
901  this->communicator_pt()->mpi_comm());
902 
903  // On very coarse meshes with larger numbers of processors, METIS
904  // occasionally returns an element_domain Vector for which a particular
905  // processor has no elements affiliated to it; the following fixes this
906 
907  // Convert element_domain to integer storage
908  Vector<int> int_element_domain(nelem);
909  for (unsigned e=0;e<nelem;e++)
910  {
911  int_element_domain[e]=element_domain[e];
912  }
913 
914  // Global storage for number of elements on each process
915  int my_number_of_elements=0;
916  Vector<int> number_of_elements(n_proc,0);
917 
918  for (unsigned e=0;e<nelem;e++)
919  {
920  if (int_element_domain[e]==rank)
921  {
922  my_number_of_elements++;
923  }
924  }
925 
926  // Communicate the correct value for each single process into
927  // the global storage vector
928  MPI_Allgather(&my_number_of_elements,1,MPI_INT,
929  &number_of_elements[0],1,MPI_INT,
930  this->communicator_pt()->mpi_comm());
931 
932  // If a process has no elements then switch an element with the
933  // process with the largest number of elements, assuming
934  // that it still has enough elements left to share
935  int max_number_of_elements=0;
936  int process_with_max_elements=0;
937  for (int d=0;d<n_proc;d++)
938  {
939  if (number_of_elements[d]==0)
940  {
941  // Find the process with maximum number of elements
942  if (max_number_of_elements<=1)
943  {
944  for (int dd=0;dd<n_proc;dd++)
945  {
946  if (number_of_elements[dd]>max_number_of_elements)
947  {
948  max_number_of_elements=number_of_elements[dd];
949  process_with_max_elements=dd;
950  }
951  }
952  }
953 
954  // Check that this number of elements is okay for sharing
955  if (max_number_of_elements<=1)
956  {
957  // Throw an error if elements can't be shared
958  std::ostringstream error_stream;
959  error_stream << "No process has more than 1 element, and\n"
960  << "at least one process has no elements!\n"
961  << "Suggest rerunning with more refinement.\n"
962  << std::endl;
963  throw OomphLibError(error_stream.str(),
964  OOMPH_CURRENT_FUNCTION,
965  OOMPH_EXCEPTION_LOCATION);
966 
967  }
968 
969  // Loop over the element domain vector and switch
970  // one value for process "process_with_max_elements" with d
971  for (unsigned e=0;e<nelem;e++)
972  {
973  if (int_element_domain[e]==process_with_max_elements)
974  {
975  int_element_domain[e]=d;
976  // Change the numbers associated with these processes
977  number_of_elements[d]++;
978  number_of_elements[process_with_max_elements]--;
979  // Reduce the number of elements available on "max" process
980  max_number_of_elements--;
981  // Inform the user that a switch has taken place
982  if (report_stats)
983  {
984  oomph_info << "INFO: Switched element domain at position " << e
985  << std::endl
986  << "from process " << process_with_max_elements
987  << " to process " << d
988  << std::endl
989  << "which was given no elements by METIS partition"
990  << std::endl;
991  }
992  // Only need to do this once for this element loop, otherwise
993  // this will take all the elements from "max" process and put them
994  // in process d, thus leaving essentially the same problem!
995  break;
996  }
997  }
998  }
999 
1000  }
1001 
1002  // Reassign new values to the element_domain vector
1003  for (unsigned e=0;e<nelem;e++)
1004  {
1005  element_domain[e]=int_element_domain[e];
1006  }
1007 
1008  unsigned count_elements=0;
1009  for (unsigned e=0; e<nelem; e++)
1010  {
1011  if(int(element_domain[e])==rank)
1012  {
1013  count_elements++;
1014  }
1015  }
1016 
1017  if (report_stats)
1018  {
1019  oomph_info << "I have " << count_elements
1020  << " elements from this partition" << std::endl << std::endl;
1021  }
1022  }
1023 
1024  //==================================================================
1025  /// (Irreversibly) prune halo(ed) elements and nodes, usually
1026  /// after another round of refinement, to get rid of
1027  /// excessively wide halo layers. Note that the current
1028  /// mesh will be now regarded as the base mesh and no unrefinement
1029  /// relative to it will be possible once this function
1030  /// has been called.
1031  //==================================================================
1033  const bool& report_stats)
1034  {
1035 
1036  // Storage for number of processors and current processor
1037  int n_proc=this->communicator_pt()->nproc();
1038 
1039  // Has the problem been distributed yet?
1041  {
1042  oomph_info
1043  << "WARNING: Problem::prune_halo_elements_and_nodes() was called on a "
1044  << "non-distributed Problem!" << std::endl;
1045  oomph_info << "Ignoring your request..." << std::endl;
1046  }
1047  else
1048  {
1049  // There are no halo layers to prune if it's a single-process job
1050  if (n_proc==1)
1051  {
1052  oomph_info << "WARNING: You've tried to prune halo layers on a problem\n"
1053  << "with only one processor: this is unnecessary.\n"
1054  << "Ignoring your request."
1055  << std::endl << std::endl;
1056  }
1057  else
1058  {
1059 
1060 #ifdef PARANOID
1061  unsigned old_ndof=ndof();
1062 #endif
1063 
1064  double t_start = 0.0;
1066  {
1067  t_start=TimingHelpers::timer();
1068  }
1069 
1070  // Call actions before distribute
1072 
1073  double t_end = 0.0;
1075  {
1076  t_end=TimingHelpers::timer();
1077  oomph_info
1078  << "Time for actions_before_distribute() in "
1079  << "Problem::prune_halo_elements_and_nodes(): "
1080  << t_end-t_start << std::endl;
1081  t_start = TimingHelpers::timer();
1082  }
1083 
1084  // Associate all elements with root in current Base mesh
1085  unsigned nel=Base_mesh_element_pt.size();
1086  std::map<GeneralisedElement*,unsigned> old_base_element_number_plus_one;
1087  std::vector<bool> old_root_is_halo_or_non_existent(nel,true);
1088  for (unsigned e=0;e<nel;e++)
1089  {
1090  // Get the base element
1092 
1093  // Does it exist locally?
1094  if (base_el_pt!=0)
1095  {
1096  // Check if it's a halo element
1097  if (!base_el_pt->is_halo())
1098  {
1099  old_root_is_halo_or_non_existent[e]=false;
1100  }
1101 
1102  // Not refineable: It's only the element iself
1103  RefineableElement* ref_el_pt=0;
1104  ref_el_pt=dynamic_cast<RefineableElement*>(base_el_pt);
1105  if (ref_el_pt==0)
1106  {
1107  old_base_element_number_plus_one[base_el_pt]=e+1;
1108  }
1109  // Refineable: Get entire tree of elements
1110  else
1111  {
1112  Vector<Tree*> tree_pt;
1113  ref_el_pt->tree_pt()->stick_all_tree_nodes_into_vector(tree_pt);
1114  unsigned ntree=tree_pt.size();
1115  for (unsigned t=0;t<ntree;t++)
1116  {
1117  old_base_element_number_plus_one[tree_pt[t]->object_pt()]=e+1;
1118  }
1119  }
1120  }
1121  }
1122 
1123 
1124 
1126  {
1127  t_end=TimingHelpers::timer();
1128  oomph_info
1129  << "Time for setup old root elements in "
1130  << "Problem::prune_halo_elements_and_nodes(): "
1131  << t_end-t_start << std::endl;
1132  t_start = TimingHelpers::timer();
1133  }
1134 
1135 
1136  // Now remember the old number of base elements
1137  unsigned nel_base_old=nel;
1138 
1139 
1140  // Prune the halo elements and nodes of the mesh(es)
1141  Vector<GeneralisedElement*> deleted_element_pt;
1142  unsigned n_mesh=nsub_mesh();
1143  if (n_mesh==0)
1144  {
1145  // Prune halo elements and nodes for the (single) global mesh
1146  mesh_pt()->prune_halo_elements_and_nodes(deleted_element_pt,
1147  doc_info,report_stats);
1148  }
1149  else
1150  {
1151  // Loop over individual submeshes and prune separately
1152  for (unsigned i_mesh=0; i_mesh<n_mesh; i_mesh++)
1153  {
1154  mesh_pt(i_mesh)->prune_halo_elements_and_nodes(deleted_element_pt,
1155  doc_info,report_stats);
1156  }
1157 
1158  // Rebuild the global mesh
1160  }
1161 
1163  {
1164  t_end=TimingHelpers::timer();
1165  oomph_info
1166  << "Total time for all mesh-level prunes in "
1167  << "Problem::prune_halo_elements_and_nodes(): "
1168  << t_end-t_start << std::endl;
1169  t_start = TimingHelpers::timer();
1170  }
1171 
1172  // Loop over all elements in newly rebuilt mesh (which contains
1173  // all element in "tree order"), find the roots
1174  // (which are either non-refineable elements or refineable elements
1175  // whose tree representations are TreeRoots)
1176  std::map<FiniteElement*,bool> root_el_done;
1177 
1178  // Vector storing vectors of pointers to new base elements associated
1179  // with the same old base element
1181  new_base_element_associated_with_old_base_element(nel_base_old);
1182 
1183  unsigned n_meshes = n_mesh;
1184  // Change the value for the number of submeshes if there is only
1185  // one mesh so that the loop below works if we have only one
1186  // mesh
1187  if (n_meshes == 0)
1188  {n_meshes = 1;}
1189 
1190  // Store which submeshes, if there are some are structured
1191  // meshes
1192  std::vector<bool> is_structured_mesh(n_meshes);
1193 
1194  // Loop over all elements in the rebuilt mesh, but make sure
1195  // that we are only looping over the structured meshes
1196  nel = 0;
1197  for (unsigned i_mesh = 0; i_mesh < n_meshes; i_mesh++)
1198  {
1199  TriangleMeshBase* tri_mesh_pt =
1200  dynamic_cast<TriangleMeshBase*>(mesh_pt(i_mesh));
1201  if (!(tri_mesh_pt!=0))
1202  {
1203  // Mark the mesh as structured mesh
1204  is_structured_mesh[i_mesh] = true;
1205  // Add the number of elements
1206  nel+=mesh_pt(i_mesh)->nelement();
1207  } // if (!(tri_mesh_pt!=0))
1208  else
1209  {
1210  // Mark the mesh as nonstructured mesh
1211  is_structured_mesh[i_mesh] = false;
1212  } // else if (!(tri_mesh_pt!=0))
1213  } // for (i_mesh < n_mesh)
1214 
1215  // Go for all the meshes (if there are submeshes)
1216  for (unsigned i_mesh = 0; i_mesh < n_meshes; i_mesh++)
1217  {
1218  // Only work with the elements in the mesh if it is an
1219  // structured mesh
1220  if (is_structured_mesh[i_mesh])
1221  {
1222  // Get the number of elements in the submesh
1223  const unsigned nele_submesh = mesh_pt(i_mesh)->nelement();
1224  for (unsigned e = 0; e < nele_submesh; e++)
1225  {
1226  // Get the element
1227  GeneralisedElement* el_pt = mesh_pt(i_mesh)->element_pt(e);
1228 
1229  // Not refineable: It's definitely a new base element
1230  RefineableElement* ref_el_pt=0;
1231  ref_el_pt=dynamic_cast<RefineableElement*>(el_pt);
1232  if (ref_el_pt==0)
1233  {
1234  unsigned old_base_el_no =
1235  old_base_element_number_plus_one[el_pt]-1;
1236  new_base_element_associated_with_old_base_element
1237  [old_base_el_no].push_back(el_pt);
1238  }
1239  // Refineable
1240  else
1241  {
1242  // Is it a tree root (after pruning)? In that case it's
1243  // a new base element
1244  if (dynamic_cast<TreeRoot*>(ref_el_pt->tree_pt()))
1245  {
1246  unsigned old_base_el_no =
1247  old_base_element_number_plus_one[el_pt]-1;
1248  new_base_element_associated_with_old_base_element
1249  [old_base_el_no].push_back(el_pt);
1250  }
1251  else
1252  {
1253  // Get associated root element
1254  FiniteElement* root_el_pt=
1255  ref_el_pt->tree_pt()->root_pt()->object_pt();
1256 
1257  if (!root_el_done[root_el_pt])
1258  {
1259  root_el_done[root_el_pt]=true;
1260  unsigned old_base_el_no=
1261  old_base_element_number_plus_one[el_pt]-1;
1262  new_base_element_associated_with_old_base_element
1263  [old_base_el_no].push_back(root_el_pt);
1264  }
1265  }
1266  }
1267  } // for (e < nele_submesh)
1268  } // if (is_structured_mesh[i_mesh])
1269  } // for (i_mesh < n_mesh)
1270 
1271  // Create a vector that stores how many new root/base elements
1272  // got spawned from each old root/base element in the global mesh
1273  Vector<unsigned> local_n_new_root(nel_base_old);
1274 #ifdef PARANOID
1275  Vector<unsigned> n_new_root_back(nel_base_old);
1276 #endif
1277  for (unsigned e=0;e<nel_base_old;e++)
1278  {
1279  local_n_new_root[e]=
1280  new_base_element_associated_with_old_base_element[e].size();
1281 
1282 #ifdef PARANOID
1283  // Backup so we can check that halo data was consistent
1284  n_new_root_back[e]=local_n_new_root[e];
1285 #endif
1286  }
1287 
1289  {
1290  t_end=TimingHelpers::timer();
1291  oomph_info
1292  << "Time for setup of new base elements in "
1293  << "Problem::prune_halo_elements_and_nodes(): "
1294  << t_end-t_start << std::endl;
1295  t_start = TimingHelpers::timer();
1296  }
1297 
1298  // Now do reduce operation to get information for all
1299  // old root/base elements -- the pruned (halo!) base elements contain
1300  // fewer associated new roots.
1301  Vector<unsigned> n_new_root(nel_base_old);
1302  MPI_Allreduce(&local_n_new_root[0],&n_new_root[0],nel_base_old,
1303  MPI_UNSIGNED,MPI_MAX,this->communicator_pt()->mpi_comm());
1304 
1305 
1307  {
1308  t_end=TimingHelpers::timer();
1309  oomph_info
1310  << "Time for allreduce in "
1311  << "Problem::prune_halo_elements_and_nodes(): "
1312  << t_end-t_start << std::endl;
1313  t_start = TimingHelpers::timer();
1314  }
1315 
1316  // Find out total number of base elements
1317  unsigned nel_base_new=0;
1318  for (unsigned e=0;e<nel_base_old;e++)
1319  {
1320  // Increment
1321  nel_base_new+=n_new_root[e];
1322 
1323 #ifdef PARANOID
1324  // If we already had data for this root previously then
1325  // the data ought to be consistent afterwards (since taking
1326  // the max of consistent numbers shouldn't change things -- this
1327  // deals with halo/haloed elements)
1328  if (!old_root_is_halo_or_non_existent[e])
1329  {
1330  if (n_new_root_back[e]!=0)
1331  {
1332  if (n_new_root_back[e]!=n_new_root[e])
1333  {
1334  std::ostringstream error_stream;
1335  error_stream
1336  << "Number of new root elements spawned from old root "
1337  << e << ": " << n_new_root[e] << "\nis not consistent"
1338  << " with previous value: " << n_new_root_back[e] << std::endl;
1339  throw OomphLibError(
1340  error_stream.str(),
1341  OOMPH_CURRENT_FUNCTION,
1342  OOMPH_EXCEPTION_LOCATION);
1343  }
1344  }
1345  }
1346 
1347 #endif
1348  }
1349 
1350  // Reset base_mesh information
1351  Base_mesh_element_pt.clear();
1352  Base_mesh_element_pt.resize(nel_base_new,0);
1354 
1355  // Now enumerate the new base/root elements consistently
1356  unsigned count=0;
1357  for (unsigned e=0;e<nel_base_old;e++)
1358  {
1359  // Old root is non-halo: Just add the new roots into the
1360  // new lookup scheme consecutively
1361  if (!old_root_is_halo_or_non_existent[e])
1362  {
1363  // Loop over new root/base element
1364  unsigned n_new_root=
1365  new_base_element_associated_with_old_base_element[e].size();
1366  for (unsigned j=0;j<n_new_root;j++)
1367  {
1368  // Store new root/base element
1369  GeneralisedElement* el_pt=
1370  new_base_element_associated_with_old_base_element[e][j];
1371  Base_mesh_element_pt[count]=el_pt;
1372  Base_mesh_element_number_plus_one[el_pt]=count+1;
1373 
1374  // Bump counter
1375  count++;
1376  }
1377  }
1378  // Old root element is halo so skip insertion (i.e. leave
1379  // entries in lookup schemes nulled) but increase counter to
1380  // ensure consistency between processors
1381  else
1382  {
1383  unsigned nskip=n_new_root[e];
1384  count+=nskip;
1385  }
1386  }
1387 
1388  // Re-setup the map between "root" element and number in global mesh
1389  // (used in the load_balance() routines)
1391 
1392 
1394  {
1395  t_end=TimingHelpers::timer();
1396  oomph_info
1397  << "Time for finishing off base mesh info "
1398  << "Problem::prune_halo_elements_and_nodes(): "
1399  << t_end-t_start << std::endl;
1400  t_start = TimingHelpers::timer();
1401  }
1402 
1403 
1404  // Call actions after distribute
1406 
1407 
1409  {
1410  t_end=TimingHelpers::timer();
1411  oomph_info
1412  << "Time for actions_after_distribute() "
1413  << "Problem::prune_halo_elements_and_nodes(): "
1414  << t_end-t_start << std::endl;
1415  t_start = TimingHelpers::timer();
1416  }
1417 
1418 
1419  // Re-assign the equation numbers (incl synchronisation if reqd)
1420 #ifdef PARANOID
1421  unsigned n_dof=assign_eqn_numbers();
1422 #else
1424 #endif
1425 
1426 
1428  {
1429  t_end=TimingHelpers::timer();
1430  oomph_info
1431  << "Time for assign_eqn_numbers() "
1432  << "Problem::prune_halo_elements_and_nodes(): "
1433  << t_end-t_start << std::endl;
1434  t_start = TimingHelpers::timer();
1435  }
1436 
1437 
1438 #ifdef PARANOID
1440  {
1441  if (n_dof!=old_ndof)
1442  {
1443  std::ostringstream error_stream;
1444  error_stream
1445  << "Number of dofs in prune_halo_elements_and_nodes() has changed "
1446  << "from " << old_ndof << " to " << n_dof << "\n"
1447  <<"Check that you've implemented any necessary actions_before/after"
1448  << "\nadapt/distribute functions, e.g. to pin redundant pressure"
1449  << " dofs etc.\n";
1450  throw OomphLibError(error_stream.str(),
1451  OOMPH_CURRENT_FUNCTION,
1452  OOMPH_EXCEPTION_LOCATION);
1453  }
1454  }
1455 #endif
1456 
1457 
1458  }
1459  }
1460 
1461  }
1462 
1463 
1464 #endif
1465 
1466 
1467 //===================================================================
1468 /// Build a single (global) mesh from a number
1469 /// of submeshes which are passed as a vector of pointers to the
1470 /// submeshes. The ordering is not necessarily optimal.
1471 //==============================================================
1473  {
1474 #ifdef PARANOID
1475  //Has a global mesh already been built
1476  if(Mesh_pt!=0)
1477  {
1478  std::string error_message =
1479  "Problem::build_global_mesh() called,\n";
1480  error_message += " but a global mesh has already been built:\n";
1481  error_message += "Problem::Mesh_pt is not zero!\n";
1482 
1483  throw OomphLibError(error_message,
1484  OOMPH_CURRENT_FUNCTION,
1485  OOMPH_EXCEPTION_LOCATION);
1486  }
1487  //Check that there are submeshes
1488  if(Sub_mesh_pt.size()==0)
1489  {
1490  std::string error_message =
1491  "Problem::build_global_mesh() called,\n";
1492  error_message += " but there are no submeshes:\n";
1493  error_message += "Problem::Sub_mesh_pt has no entries\n";
1494 
1495  throw OomphLibError(error_message,
1496  OOMPH_CURRENT_FUNCTION,
1497  OOMPH_EXCEPTION_LOCATION);
1498  }
1499 #endif
1500 
1501  //Create an empty mesh
1502  Mesh_pt = new Mesh();
1503 
1504  //Call the rebuild function to construct the mesh
1506  }
1507 
1508 //====================================================================
1509 /// If one of the submeshes has changed (e.g. by
1510 /// mesh adaptation) we need to update the global mesh.
1511 /// \b Note: The nodes boundary information refers to the
1512 /// boundary numbers within the submesh!
1513 /// N.B. This is essentially the same function as the Mesh constructor
1514 /// that assembles a single global mesh from submeshes
1515 //=====================================================================
1517  {
1518  // Use the function in mesh to merge the submeshes into this one
1520  }
1521 
1522 
1523 //================================================================
1524 /// Add a timestepper to the problem. The function will automatically
1525 /// create or resize the Time object so that it contains the appropriate
1526 /// number of levels of storage.
1527 //================================================================
1528  void Problem::add_time_stepper_pt(TimeStepper* const &time_stepper_pt)
1529  {
1530  //Add the timestepper to the vector
1531  Time_stepper_pt.push_back(time_stepper_pt);
1532 
1533  //Find the number of timesteps required by the timestepper
1534  unsigned ndt = time_stepper_pt->ndt();
1535 
1536  //If time has not been allocated, create time object with the
1537  //required number of time steps
1538  if(Time_pt==0)
1539  {
1540  Time_pt = new Time(ndt);
1541  oomph_info << "Created Time with " << ndt << " timesteps" << std::endl;
1542  }
1543  else
1544  {
1545  //If the required number of time steps is greater than currently stored
1546  //resize the time storage
1547  if(ndt > Time_pt->ndt())
1548  {
1549  Time_pt->resize(ndt);
1550  oomph_info << "Resized Time to include " << ndt << " timesteps"
1551  << std::endl;
1552  }
1553  //Otherwise report that we are OK
1554  else
1555  {
1556  oomph_info << "Time object already has storage for " << ndt
1557  << " timesteps" << std::endl;
1558  }
1559  }
1560 
1561  //Pass the pointer to time to the timestepper
1562  time_stepper_pt->time_pt() = Time_pt;
1563  }
1564 
1565 //================================================================
1566 /// Set the explicit time stepper for the problem and also
1567 /// ensure that a time object has been created.
1568 //================================================================
1570  &explicit_time_stepper_pt)
1571  {
1572  //Set the explicit time stepper
1574 
1575  //If time has not been allocated, create time object with the
1576  //required number of time steps
1577  if(Time_pt==0)
1578  {
1579  Time_pt = new Time(0);
1580  oomph_info << "Created Time with storage for no previous timestep"
1581  << std::endl;
1582  }
1583  else
1584  {
1585  oomph_info << "Time object already exists " << std::endl;
1586  }
1587 
1588  }
1589 
1590 
1591 #ifdef OOMPH_HAS_MPI
1592 
1593 //================================================================
1594 /// Set default first and last elements for parallel assembly
1595 /// of non-distributed problem.
1596 //================================================================
1598  {
1600  {
1601 
1602  // Minimum number of elements per processor if there are fewer elements
1603  // than processors
1604  unsigned min_el=10;
1605 
1606  // Resize and make default assignments
1607  int n_proc=this->communicator_pt()->nproc();
1608  unsigned n_elements=Mesh_pt->nelement();
1609  First_el_for_assembly.resize(n_proc,0);
1610  Last_el_plus_one_for_assembly.resize(n_proc,0);
1611 
1612  // In the absence of any better knowledge distribute work evenly
1613  // over elements
1614  unsigned range=0;
1615  unsigned lo_proc=0;
1616  unsigned hi_proc=n_proc-1;
1617  if (int(n_elements)>=n_proc)
1618  {
1619  range=unsigned(double(n_elements)/double(n_proc));
1620  }
1621  else
1622  {
1623  range=min_el;
1624  lo_proc=0;
1625  hi_proc=unsigned(double(n_elements)/double(min_el));
1626  }
1627 
1628  for (int p=lo_proc;p<=int(hi_proc);p++)
1629  {
1630  First_el_for_assembly[p] = p*range;
1631 
1632  unsigned last_el_plus_one=(p+1)*range;
1633  if (last_el_plus_one>n_elements) last_el_plus_one=n_elements;
1634  Last_el_plus_one_for_assembly[p] = last_el_plus_one;
1635  }
1636 
1637  // Last one needs to incorporate any dangling elements
1638  if (int(n_elements)>=n_proc)
1639  {
1640  Last_el_plus_one_for_assembly[n_proc-1] = n_elements;
1641  }
1642 
1643  // Doc
1644  if (n_proc>1)
1645  {
1646 
1648  {
1649  oomph_info << "Problem is not distributed. Parallel assembly of "
1650  << "Jacobian uses default partitioning: "<< std::endl;
1651  for (int p=0;p<n_proc;p++)
1652  {
1654  {
1655  oomph_info << "Proc " << p << " assembles from element "
1656  << First_el_for_assembly[p] << " to "
1657  << Last_el_plus_one_for_assembly[p]-1 << " \n";
1658  }
1659  else
1660  {
1661  oomph_info << "Proc " << p << " assembles no elements\n";
1662  }
1663  }
1664  }
1665  }
1666  }
1667 
1668  }
1669 
1670 
1671 //=======================================================================
1672 /// Helper function to re-assign the first and last elements to be
1673 /// assembled by each processor during parallel assembly for
1674 /// non-distributed problem.
1675 //=======================================================================
1677  {
1678  // Wait until all processes have completed/timed their assembly
1679  MPI_Barrier(this->communicator_pt()->mpi_comm());
1680 
1681  // Storage for number of processors and current processor
1682  int n_proc=this->communicator_pt()->nproc();
1683  int rank=this->communicator_pt()->my_rank();
1684 
1685  // Don't bother to update if we've got fewer elements than
1686  // processors
1687  unsigned nel=Elemental_assembly_time.size();
1688  if (int(nel)<n_proc)
1689  {
1690  oomph_info << "Not re-computing distribution of elemental assembly\n"
1691  << "because there are fewer elements than processors\n";
1692  return;
1693  }
1694 
1695  // Setup vectors storing the number of element timings to be sent
1696  // and the offset in the final vector
1697  Vector<int> receive_count(n_proc);
1698  Vector<int> displacement(n_proc);
1699  int offset=0;
1700  for (int p=0;p<n_proc;p++)
1701  {
1702  // Default distribution of labour
1703  unsigned el_lo = First_el_for_assembly[p];
1704  unsigned el_hi = Last_el_plus_one_for_assembly[p]-1;
1705 
1706  // Number of timings to be sent and offset from start in
1707  // final vector
1708  receive_count[p]=el_hi-el_lo+1;
1709  displacement[p]=offset;
1710  offset+=el_hi-el_lo+1;
1711  }
1712 
1713  // Make temporary c-style array to avoid over-writing in Gatherv below
1714  double* el_ass_time = new double[nel];
1715  for (unsigned e=0;e<nel;e++)
1716  {
1717  el_ass_time[e]=Elemental_assembly_time[e];
1718  }
1719 
1720  // Gather timings on root processor
1721  unsigned nel_local =Last_el_plus_one_for_assembly[rank]-1
1722  -First_el_for_assembly[rank]+1;
1723  MPI_Gatherv(
1724  &el_ass_time[First_el_for_assembly[rank]],nel_local,MPI_DOUBLE,
1725  &Elemental_assembly_time[0],&receive_count[0],&displacement[0],
1726  MPI_DOUBLE,0,this->communicator_pt()->mpi_comm());
1727  delete[] el_ass_time;
1728 
1729  // Vector of first and last elements for each processor
1730  Vector<Vector <int> > first_and_last_element(n_proc);
1731  for(int p=0;p<n_proc;p++) {first_and_last_element[p].resize(2);}
1732 
1733  // Re-distribute work
1734  if (rank==0)
1735  {
1737  {
1738  oomph_info
1739  << std::endl
1740  << "Re-assigning distribution of element assembly over processors:"
1741  << std::endl;
1742  }
1743 
1744  // Get total assembly time
1745  double total=0.0;
1746  unsigned n_elements=Mesh_pt->nelement();
1747  for (unsigned e=0;e<n_elements;e++)
1748  {
1749  total+=Elemental_assembly_time[e];
1750  }
1751 
1752  // Target load per processor
1753  double target_load=total/double(n_proc);
1754 
1755  // We're on the root processor: Always start with the first element
1756  int proc=0;
1757  first_and_last_element[0][0] = 0;
1758 
1759  // Highest element we can help ourselves to if we want to leave
1760  // at least one element for all subsequent processors
1761  unsigned max_el_avail=n_elements-n_proc;
1762 
1763  // Initialise total work allocated
1764  total=0.0;
1765  for (unsigned e=0;e<n_elements;e++)
1766  {
1767  total+=Elemental_assembly_time[e];
1768 
1769  //Once we have reached the target load or we've used up practically
1770  // all the elements...
1771  if ((total>target_load)||(e==max_el_avail))
1772  {
1773 
1774  // Last element for current processor
1775  first_and_last_element[proc][1]=e;
1776 
1777  //Provided that we are not on the last processor
1778  if (proc<(n_proc-1))
1779  {
1780  // Set first element for next one
1781  first_and_last_element[proc+1][0]=e+1;
1782 
1783  // Move on to the next processor
1784  proc++;
1785  }
1786 
1787  // Can have one more...
1788  max_el_avail++;
1789 
1790  // Re-initialise the time
1791  total=0.0;
1792  } // end of test for "total exceeds target"
1793  }
1794 
1795 
1796  // Last element for last processor
1797  first_and_last_element[n_proc-1][1]=n_elements-1;
1798 
1799 
1800  // The following block should probably be paranoidified away
1801  // but we've screwed the logic up so many times that I feel
1802  // it's safer to keep it...
1803  bool wrong=false;
1804  std::ostringstream error_stream;
1805  for (int p=0;p<n_proc-1;p++)
1806  {
1807  unsigned first_of_current=first_and_last_element[p][0];
1808  unsigned last_of_current=first_and_last_element[p][1];
1809  if (first_of_current>last_of_current)
1810  {
1811  wrong=true;
1812  error_stream << "Error: First/last element of proc " << p << ": "
1813  << first_of_current << " " << last_of_current
1814  << std::endl;
1815  }
1816  unsigned first_of_next=first_and_last_element[p+1][0];
1817  if (first_of_next!=(last_of_current+1))
1818  {
1819  wrong=true;
1820  error_stream << "Error: First element of proc " << p+1 << ": "
1821  << first_of_next << " and last element of proc "
1822  << p << ": " << last_of_current
1823  << std::endl;
1824  }
1825  }
1826  if (wrong)
1827  {
1828  throw OomphLibError(error_stream.str(),
1829  OOMPH_CURRENT_FUNCTION,
1830  OOMPH_EXCEPTION_LOCATION);
1831  }
1832 
1833 
1834  // THIS TIDY UP SHOULD NO LONGER BE REQUIRED AND CAN GO AT SOME POINT
1835 
1836  // //If we haven't got to the end of the processor list then
1837  // //need to shift things about slightly because the processors
1838  // //at the end will be empty.
1839  // //This can occur when you have very fast assembly times and the
1840  // //rounding errors mean that the targets are achieved before all processors
1841  // //have been visited.
1842  // //Happens a lot when you massively oversubscribe the CPUs (which was
1843  // //only ever for testing!)
1844  // if (proc!=n_proc-1)
1845  // {
1846  // oomph_info
1847  // << "First pass did not allocate elements on every processor\n";
1848  // oomph_info <<
1849  // "Moving elements so that each processor has at least one\n";
1850 
1851  // //Work out number of empty processos
1852  // unsigned n_empty_processors = n_proc - proc + 1;
1853 
1854  // //Loop over the processors that do have elements
1855  // //and work out how many we need to steal elements from
1856  // unsigned n_element_on_processors=0;
1857  // do
1858  // {
1859  // //Step down the processors
1860  // --proc;
1861  // //Add the current processor to the number of empty processors
1862  // //because the elements have to be shared between processors
1863  // //including the one(s) on which they are currently stored.
1864  // ++n_empty_processors;
1865  // n_element_on_processors +=
1866  // (first_and_last_element[proc][1] -
1867  // first_and_last_element[proc][0] + 1);
1868  // }
1869  // while(n_element_on_processors < n_empty_processors);
1870 
1871  // //Should now be able to put one element on each processor
1872  // //Start from the end and do so
1873  // unsigned current_element = n_elements-1;
1874  // for(int p=n_proc-1;p>proc;p--)
1875  // {
1876  // first_and_last_element[p][1] = current_element;
1877  // first_and_last_element[p][0] = --current_element;
1878  // }
1879 
1880  // //Now for the last processor we touched, just adjust the final value
1881  // first_and_last_element[proc][1] = current_element;
1882  // }
1883  // //Otherwise just put the rest of the elements on the final
1884  // //processor
1885  // else
1886  // {
1887  // // Last one
1888  // first_and_last_element[n_proc-1][1]=n_elements-1;
1889  // }
1890 
1891 
1892  // END PRESUMED-TO-BE-UNNECESSARY BLOCK...
1893 
1894 
1895 
1896  //Now communicate the information
1897 
1898  //Set local informationt for this (root) processor
1899  First_el_for_assembly[0]= first_and_last_element[0][0];
1900  Last_el_plus_one_for_assembly[0] = first_and_last_element[0][1] + 1;
1901 
1903  {
1904  oomph_info
1905  << "Processor " << 0 << " assembles Jacobians"
1906  << " from elements " << first_and_last_element[0][0] << " to "
1907  << first_and_last_element[0][1] << " "
1908  << std::endl;
1909  }
1910 
1911  //Only now can we send the information to the other processors
1912  for(int p=1;p<n_proc;++p)
1913  {
1914 
1915  MPI_Send(&first_and_last_element[p][0],2,MPI_INT,p,
1916  0,this->communicator_pt()->mpi_comm());
1917 
1918 
1920  {
1921  oomph_info
1922  << "Processor " << p << " assembles Jacobians"
1923  << " from elements " << first_and_last_element[p][0] << " to "
1924  << first_and_last_element[p][1] << " "
1925  << std::endl;
1926  }
1927  }
1928  }
1929  // Receive first and last element from root on non-master processors
1930  else
1931  {
1932  Vector<int> aux(2);
1933  MPI_Status status;
1934  MPI_Recv(&aux[0],2,MPI_INT,0,0,
1935  this->communicator_pt()->mpi_comm(),&status);
1936  First_el_for_assembly[rank]=aux[0];
1937  Last_el_plus_one_for_assembly[rank]=aux[1]+1;
1938  }
1939 
1940  // Wipe all others
1941  for (int p=0;p<n_proc;p++)
1942  {
1943  if (p!=rank)
1944  {
1945  First_el_for_assembly[p]=0;
1947  }
1948  }
1949 
1950  // The equations assembled by this processor may have changed so
1951  // we must resize the sparse assemble with arrays previous allocation
1953  }
1954 
1955 #endif
1956 
1957 //================================================================
1958 /// Assign all equation numbers for problem: Deals with global
1959 /// data (= data that isn't attached to any elements) and then
1960 /// does the equation numbering for the elements. Bool argument
1961 /// can be set to false to ignore assigning local equation numbers
1962 /// (necessary in the parallel implementation of locate_zeta
1963 /// between multiple meshes).
1964 //================================================================
1965  unsigned long Problem::assign_eqn_numbers(const bool& assign_local_eqn_numbers)
1966  {
1967  //Check that the global mesh has been build
1968 #ifdef PARANOID
1969  if(Mesh_pt==0)
1970  {
1971  std::ostringstream error_stream;
1972  error_stream <<
1973  "Global mesh does not exist, so equation numbers cannot be assigned.\n";
1974  //Check for sub meshes
1975  if(nsub_mesh()==0)
1976  {
1977  error_stream << "There aren't even any sub-meshes in the Problem.\n"
1978  << "You can set the global mesh directly by using\n"
1979  << "Problem::mesh_pt() = my_mesh_pt;\n"
1980  << "OR you can use Problem::add_sub_mesh(mesh_pt); "
1981  << "to add a sub mesh.\n";
1982  }
1983  else
1984  {
1985  error_stream << "There are " << nsub_mesh() << " sub-meshes.\n";
1986  }
1987  error_stream <<
1988  "You need to call Problem::build_global_mesh() to create a global mesh\n"
1989  << "from the sub-meshes.\n\n";
1990 
1991  throw OomphLibError(error_stream.str(),
1992  OOMPH_CURRENT_FUNCTION,
1993  OOMPH_EXCEPTION_LOCATION);
1994  }
1995 #endif
1996 
1997  // Number of submeshes
1998  unsigned n_sub_mesh=Sub_mesh_pt.size();
1999 
2000 #ifdef OOMPH_HAS_MPI
2001 
2002  // Storage for number of processors
2003  int n_proc=this->communicator_pt()->nproc();
2004 
2005 
2006  if (n_proc>1)
2007  {
2008  // Force re-analysis of time spent on assembly each
2009  // elemental Jacobian
2011  Elemental_assembly_time.clear();
2012  }
2013  else
2014  {
2016  }
2017 
2018  // Re-distribution of elements over processors during assembly
2019  // must be recomputed
2021  {
2022  // Set default first and last elements for parallel assembly
2023  // of non-distributed problem.
2025  }
2026 
2027 #endif
2028 
2029 
2030  double t_start = 0.0;
2032  {
2033  t_start=TimingHelpers::timer();
2034  }
2035 
2036  // Loop over all elements in the mesh and set up any additional
2037  // dependencies that they may have (e.g. storing the geometric
2038  // Data, i.e. Data that affects an element's shape in elements
2039  // with algebraic node-update functions
2040  unsigned nel=Mesh_pt->nelement();
2041  for (unsigned e=0;e<nel;e++)
2042  {
2044  }
2045 
2046 #ifdef OOMPH_HAS_MPI
2047  // Complete setup of dependencies for external halo elements too
2048  unsigned n_mesh=this->nsub_mesh();
2049  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
2050  {
2051  for (int iproc=0;iproc<n_proc;iproc++)
2052  {
2053  unsigned n_ext_halo_el=mesh_pt(i_mesh)->nexternal_halo_element(iproc);
2054  for (unsigned e=0;e<n_ext_halo_el;e++)
2055  {
2056  mesh_pt(i_mesh)->external_halo_element_pt(iproc,e)
2058  }
2059  }
2060  }
2061 #endif
2062 
2063 
2064 
2065 
2066  double t_end = 0.0;
2068  {
2069  t_end = TimingHelpers::timer();
2070  oomph_info
2071  << "Time for complete setup of dependencies in assign_eqn_numbers: "
2072  << t_end-t_start << std::endl;
2073  }
2074 
2075 
2076  // Initialise number of dofs for reserve below
2077  unsigned n_dof=0;
2078 
2079  // Potentially loop over remainder of routine, possible re-visiting all those
2080  // parts that must be redone, following the removal of duplicate
2081  // external halo data.
2082  for (unsigned loop_count=0;loop_count<2;loop_count++)
2083  {
2084  //(Re)-set the dof pointer to zero length because entries are
2085  //pushed back onto it -- if it's not reset here then we get into
2086  //trouble during mesh refinement when we reassign all dofs
2087  Dof_pt.resize(0);
2088 
2089  // Reserve from previous allocation if we're going around again
2090  Dof_pt.reserve(n_dof);
2091 
2092  //Reset the equation number
2093  unsigned long equation_number=0;
2094 
2095  //Now set equation numbers for the global Data
2096  unsigned Nglobal_data = nglobal_data();
2097  for(unsigned i=0;i<Nglobal_data;i++)
2098  {Global_data_pt[i]->assign_eqn_numbers(equation_number,Dof_pt);}
2099 
2101  {
2102  t_start = TimingHelpers::timer();
2103  }
2104 
2105  //Call assign equation numbers on the global mesh
2107 
2108  // Deal with the spine meshes additional numbering
2109  //If there is only one mesh
2110  if(n_sub_mesh==0)
2111  {
2112  if(SpineMesh* const spine_mesh_pt = dynamic_cast<SpineMesh*>(Mesh_pt))
2113  {
2114  n_dof = spine_mesh_pt->assign_global_spine_eqn_numbers(Dof_pt);
2115  }
2116  }
2117  //Otherwise loop over the sub meshes
2118  else
2119  {
2120  //Assign global equation numbers first
2121  for(unsigned i=0;i<n_sub_mesh;i++)
2122  {
2123  if(SpineMesh* const spine_mesh_pt =
2124  dynamic_cast<SpineMesh*>(Sub_mesh_pt[i]))
2125  {
2126  n_dof = spine_mesh_pt->assign_global_spine_eqn_numbers(Dof_pt);
2127  }
2128  }
2129  }
2130 
2132  {
2133  t_end = TimingHelpers::timer();
2134  oomph_info
2135  << "Time for assign_global_eqn_numbers in assign_eqn_numbers: "
2136  << t_end-t_start << std::endl;
2137  t_start = TimingHelpers::timer();
2138  }
2139 
2140 
2141 #ifdef OOMPH_HAS_MPI
2142 
2143  // reset previous allocation
2145 
2146  // Only synchronise if the problem has actually been
2147  // distributed.
2149  {
2150  // Synchronise the equation numbers and return the total
2151  // number of degrees of freedom in the overall problem
2152  // Do not assign local equation numbers -- we're doing this
2153  // below.
2154  n_dof=synchronise_eqn_numbers(false);
2155  }
2156  // ..else just setup the Dof_distribution_pt
2157  // NOTE: this is setup by synchronise_eqn_numbers(...)
2158  // if Problem_has_been_distributed
2159  else
2160 #endif
2161  {
2163  }
2164 
2166  {
2167  t_end = TimingHelpers::timer();
2168  oomph_info
2169  << "Time for Problem::synchronise_eqn_numbers in "
2170  << "Problem::assign_eqn_numbers: "
2171  << t_end-t_start << std::endl;
2172  }
2173 
2174 
2175 #ifdef OOMPH_HAS_MPI
2176 
2177 
2178  // Now remove duplicate data in external halo elements
2180  {
2181 
2183  {
2184  t_start = TimingHelpers::timer();
2185  }
2186 
2187  // Monitor if we've actually changed anything
2188  bool actually_removed_some_data=false;
2189 
2190  // Only do it once!
2191  if (loop_count==0)
2192  {
2193  if (n_sub_mesh==0)
2194  {
2195  remove_duplicate_data(Mesh_pt,actually_removed_some_data);
2196  }
2197  else
2198  {
2199  for (unsigned i=0;i<n_sub_mesh;i++)
2200  {
2201  bool tmp_actually_removed_some_data=false;
2203  tmp_actually_removed_some_data);
2204  if (tmp_actually_removed_some_data) actually_removed_some_data=true;
2205  }
2206  }
2207  }
2208 
2209 
2211  {
2212  t_end = TimingHelpers::timer();
2213  std::stringstream tmp;
2214  tmp << "Time for calls to Problem::remove_duplicate_data in "
2215  << "Problem::assign_eqn_numbers: "
2216  << t_end-t_start << " ; have ";
2217  if (!actually_removed_some_data)
2218  {
2219  tmp << " not ";
2220  }
2221  tmp << " removed some/any data.\n";
2222  oomph_info << tmp.str();
2223  t_start=TimingHelpers::timer();
2224  }
2225 
2226  // Break out of the loop if we haven't done anything here.
2227  unsigned status=0;
2228  if (actually_removed_some_data) status=1;
2229 
2230  // Allreduce to check if anyone has removed any data
2231  unsigned overall_status=0;
2232  MPI_Allreduce(&status,&overall_status,1,
2233  MPI_UNSIGNED,MPI_MAX,this->communicator_pt()->mpi_comm());
2234 
2235 
2237  {
2238  t_end = TimingHelpers::timer();
2239  std::stringstream tmp;
2240  tmp << "Time for MPI_Allreduce after Problem::remove_duplicate_data in "
2241  << "Problem::assign_eqn_numbers: "
2242  << t_end-t_start << std::endl;
2243  oomph_info << tmp.str();
2244  t_start=TimingHelpers::timer();
2245  }
2246 
2247  // Bail out if we haven't done anything here
2248  if (overall_status!=1)
2249  {
2250  break;
2251  }
2252 
2253  // Big tidy up: Remove null pointers from halo/haloed node storage
2254  // for all meshes (this involves comms and therefore must be
2255  // performed outside loop over meshes so the all-to-all is only
2256  // done once)
2258 
2259  // Time it...
2261  {
2262  double t_end = TimingHelpers::timer();
2263  oomph_info
2264  << "Total time for "
2265  << "Problem::remove_null_pointers_from_external_halo_node_storage(): "
2266  << t_end-t_start << std::endl;
2267  }
2268 
2269  }
2270  else
2271  {
2272  // Problem not distributed; no need for another loop
2273  break;
2274  }
2275 
2276 #else
2277 
2278  // Serial run: Again no need for a second loop
2279  break;
2280 
2281 #endif
2282 
2283  } // end of loop over fcts that need to be re-executed if
2284  // we've removed duplicate external data
2285 
2286 
2287  // Resize the sparse assemble with arrays previous allocation
2289 
2290 
2292  {
2293  t_start = TimingHelpers::timer();
2294  }
2295 
2296  // Finally assign local equations
2297  if (assign_local_eqn_numbers)
2298  {
2299  if (n_sub_mesh==0)
2300  {
2302  }
2303  else
2304  {
2305  for (unsigned i=0;i<n_sub_mesh;i++)
2306  {
2307  Sub_mesh_pt[i]->
2308  assign_local_eqn_numbers(Store_local_dof_pt_in_elements);
2309  }
2310  }
2311  }
2312 
2314  {
2315  t_end = TimingHelpers::timer();
2316  oomph_info
2317  << "Total time for all Mesh::assign_local_eqn_numbers in "
2318  << "Problem::assign_eqn_numbers: "
2319  << t_end-t_start << std::endl;
2320  }
2321 
2322 
2323  // and return the total number of DOFs
2324  return n_dof;
2325 
2326  }
2327 //================================================================
2328 /// \short Function to describe the dofs in terms of the global
2329 /// equation number, i.e. what type of value (nodal value of
2330 /// a Node; value in a Data object; value of internal Data in an
2331 /// element; etc) is the unknown with a certain global equation number.
2332 /// Output stream defaults to oomph_info.
2333 //================================================================
2334  void Problem::describe_dofs(std::ostream& out) const
2335  {
2336  //Check that the global mesh has been build
2337 #ifdef PARANOID
2338  if(Mesh_pt==0)
2339  {
2340  std::ostringstream error_stream;
2341  error_stream <<
2342  "Global mesh does not exist, so equation numbers cannot be found.\n";
2343  //Check for sub meshes
2344  if(nsub_mesh()==0)
2345  {
2346  error_stream << "There aren't even any sub-meshes in the Problem.\n"
2347  << "You can set the global mesh directly by using\n"
2348  << "Problem::mesh_pt() = my_mesh_pt;\n"
2349  << "OR you can use Problem::add_sub_mesh(mesh_pt); "
2350  << "to add a sub mesh.\n";
2351  }
2352  else
2353  {
2354  error_stream << "There are " << nsub_mesh() << " sub-meshes.\n";
2355  }
2356  error_stream <<
2357  "You need to call Problem::build_global_mesh() to create a global mesh\n"
2358  << "from the sub-meshes.\n\n";
2359 
2360  throw OomphLibError(error_stream.str(),
2361  OOMPH_CURRENT_FUNCTION,
2362  OOMPH_EXCEPTION_LOCATION);
2363  }
2364 #endif
2365 
2366  out << "Although this program will describe the degrees of freedom in the \n"
2367  << "problem, it will do so using the typedef for the elements. This is \n"
2368  << "not neccesarily human readable, but there is a solution.\n"
2369  << "Pipe your program's output through c++filt, with the argument -t.\n"
2370  << "e.g. \"./two_d_multi_poisson | c++filt -t > ReadableOutput.txt\".\n "
2371  << "(Disregarding the quotes)\n\n\n";
2372 
2373  out << "Classifying Global Equation Numbers" << std::endl;
2374  out << std::string(80,'-') << std::endl;
2375 
2376  // Number of submeshes
2377  unsigned n_sub_mesh=Sub_mesh_pt.size();
2378 
2379  // Classify Global dofs
2380  unsigned Nglobal_data = nglobal_data();
2381  for(unsigned i=0;i<Nglobal_data;i++)
2382  {
2383  std::stringstream conversion;
2384  conversion <<" in Global Data "<<i<<".";
2385  std::string in(conversion.str());
2386  Global_data_pt[i]->describe_dofs(out,in);
2387  }
2388 
2389  // Put string in limiting scope.
2390  {
2391  // Descend into assignment for mesh.
2392  std::string in(" in Problem's Only Mesh.");
2393  Mesh_pt->describe_dofs(out,in);
2394  }
2395 
2396  // Deal with the spine meshes additional numbering:
2397  // If there is only one mesh:
2398  if(n_sub_mesh==0)
2399  {
2400  if(SpineMesh* const spine_mesh_pt = dynamic_cast<SpineMesh*>(Mesh_pt))
2401  {
2402  std::string in(" in Problem's Only SpineMesh.");
2403  spine_mesh_pt->describe_spine_dofs(out,in);
2404  }
2405  }
2406  //Otherwise loop over the sub meshes
2407  else
2408  {
2409  //Assign global equation numbers first
2410  for(unsigned i=0;i<n_sub_mesh;i++)
2411  {
2412  if(SpineMesh* const spine_mesh_pt =
2413  dynamic_cast<SpineMesh*>(Sub_mesh_pt[i]))
2414  {
2415  std::stringstream conversion;
2416  conversion <<" in Sub-SpineMesh "<<i<<".";
2417  std::string in(conversion.str());
2418  spine_mesh_pt->describe_spine_dofs(out,in);
2419  }// end if.
2420  }// end for.
2421  }// end else.
2422 
2423 
2424  out << std::string(80,'\\') << std::endl;
2425  out << std::string(80,'\\') << std::endl;
2426  out << std::string(80,'\\') << std::endl;
2427  out << "Classifying global eqn numbers in terms of elements." << std::endl;
2428  out << std::string(80,'-') << std::endl;
2429  out << "Eqns | Source" << std::endl;
2430  out << std::string(80,'-') << std::endl;
2431 
2432  if (n_sub_mesh==0)
2433  {
2434  std::string in(" in Problem's Only Mesh.");
2435  Mesh_pt->describe_local_dofs(out,in);
2436  }
2437  else
2438  {
2439  for (unsigned i=0;i<n_sub_mesh;i++)
2440  {
2441  std::stringstream conversion;
2442  conversion <<" in Sub-Mesh "<<i<<".";
2443  std::string in(conversion.str());
2444  Sub_mesh_pt[i]->describe_local_dofs(out,in);
2445  }// End for
2446  }// End else
2447  } // End problem::describe_dofs(...)
2448 
2449 
2450 //================================================================
2451 /// Get the vector of dofs, i.e. a vector containing the current
2452 /// values of all unknowns.
2453 //================================================================
2455  {
2456  //Find number of dofs
2457  const unsigned long n_dof = ndof();
2458 
2459  //Resize the vector
2460  dofs.build(Dof_distribution_pt,0.0);
2461 
2462  //Copy dofs into vector
2463  for(unsigned long l=0;l<n_dof;l++)
2464  {
2465  dofs[l] = *Dof_pt[l];
2466  }
2467  }
2468 
2469  /// Get history values of dofs in a double vector.
2470  void Problem::get_dofs(const unsigned& t, DoubleVector& dofs) const
2471  {
2472 #ifdef PARANOID
2473  if(distributed())
2474  {
2475  throw OomphLibError("Not designed for distributed problems",
2476  OOMPH_EXCEPTION_LOCATION,
2477  OOMPH_CURRENT_FUNCTION);
2478  // might work, not sure
2479  }
2480 #endif
2481 
2482  // Resize the vector
2483  dofs.build(Dof_distribution_pt, 0.0);
2484 
2485  // First deal with global data
2486  unsigned Nglobal_data = nglobal_data();
2487  for(unsigned i=0; i<Nglobal_data; i++)
2488  {
2489  for(unsigned j=0, nj=Global_data_pt[i]->nvalue(); j<nj; j++)
2490  {
2491  // For each data get the equation number and copy out the value.
2492  int eqn_number = Global_data_pt[i]->eqn_number(j);
2493  if(eqn_number >= 0)
2494  {
2495  dofs[eqn_number] = Global_data_pt[i]->value(t, j);
2496  }
2497  }
2498  }
2499 
2500  // Next element internal data
2501  for(unsigned i=0, ni=mesh_pt()->nelement(); i<ni; i++)
2502  {
2503 
2504  GeneralisedElement* ele_pt = mesh_pt()->element_pt(i);
2505  for(unsigned j=0, nj=ele_pt->ninternal_data(); j<nj; j++)
2506  {
2507 
2508  Data* d_pt = ele_pt->internal_data_pt(j);
2509  for(unsigned k=0, nk=d_pt->nvalue(); k<nk; k++)
2510  {
2511 
2512  int eqn_number = d_pt->eqn_number(k);
2513  if(eqn_number >= 0)
2514  {
2515  dofs[eqn_number] = d_pt->value(t, k);
2516  }
2517  }
2518  }
2519  }
2520 
2521  // Now the nodes
2522  for(unsigned i=0, ni=mesh_pt()->nnode(); i<ni; i++)
2523  {
2524  Node* node_pt = mesh_pt()->node_pt(i);
2525  for(unsigned j=0, nj=node_pt->nvalue(); j<nj; j++)
2526  {
2527  // For each node get the equation number and copy out the value.
2528  int eqn_number = node_pt->eqn_number(j);
2529  if(eqn_number >= 0)
2530  {
2531  dofs[eqn_number] = node_pt->value(t, j);
2532 
2533  }
2534  }
2535  }
2536  }
2537 
2538 
2539 #ifdef OOMPH_HAS_MPI
2540 
2541 //=======================================================================
2542 /// Private helper function to remove repeated data
2543 /// in external haloed elements associated with specified mesh.
2544 /// Bool is true if some data was removed -- this usually requires
2545 /// re-running through certain parts of the equation numbering procedure.
2546 //======================================================================
2547  void Problem::remove_duplicate_data(Mesh* const &mesh_pt,
2548  bool& actually_removed_some_data)
2549  {
2550 
2551 // // // Taken out again by MH -- clutters up output
2552 // // Doc timings if required
2553 // double t_start=0.0;
2554 // if (Global_timings::Doc_comprehensive_timings)
2555 // {
2556 // t_start=TimingHelpers::timer();
2557 // }
2558 
2559  int n_proc=this->communicator_pt()->nproc();
2560  int my_rank=this->communicator_pt()->my_rank();
2561 
2562  // Initialise
2563  actually_removed_some_data=false;
2564 
2565  // Each individual container of external halo nodes has unique
2566  // nodes/equation numbers, but there may be some duplication between
2567  // two or more different containers; the following code checks for this
2568  // and removes the duplication by overwriting any data point with an already
2569  // existing eqn number with the original data point which had the eqn no.
2570 
2571  // // Storage for existing nodes, enumerated by first non-negative
2572  // // global equation number
2573  // unsigned n_dof=ndof();
2574 
2575  // Note: This used to be
2576  // Vector<Node*> global_node_pt(n_dof,0);
2577  // but this is a total killer! Memory allocation is extremely
2578  // costly and only relatively few entries are used so use
2579  // map:
2580  std::map<unsigned,Node*> global_node_pt;
2581 
2582  // Only do each retained node once
2583  std::map<Node*,bool> node_done;
2584 
2585  // Loop over existing "normal" elements in mesh
2586  unsigned n_element=mesh_pt->nelement();
2587  for (unsigned e=0;e<n_element;e++)
2588  {
2589  FiniteElement* el_pt = dynamic_cast<FiniteElement*>(
2590  mesh_pt->element_pt(e));
2591  if (el_pt!=0)
2592  {
2593  // Loop over nodes
2594  unsigned n_node=el_pt->nnode();
2595  for (unsigned j=0;j<n_node;j++)
2596  {
2597  Node* nod_pt=el_pt->node_pt(j);
2598 
2599  // Have we already done the node?
2600  if (!node_done[nod_pt])
2601  {
2602  node_done[nod_pt]=true;
2603 
2604  // Loop over values stored at node (if any) to find
2605  // the first non-negative eqn number
2606  unsigned first_non_negative_eqn_number_plus_one=0;
2607  unsigned n_val=nod_pt->nvalue();
2608  for (unsigned i_val=0;i_val<n_val;i_val++)
2609  {
2610  int eqn_no=nod_pt->eqn_number(i_val);
2611  if (eqn_no>=0)
2612  {
2613  first_non_negative_eqn_number_plus_one=eqn_no+1;
2614  break;
2615  }
2616  }
2617 
2618  // If we haven't found a non-negative eqn number check
2619  // eqn numbers associated with solid data (if any)
2620  if (first_non_negative_eqn_number_plus_one==0)
2621  {
2622  // Is it a solid node?
2623  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
2624  if (solid_nod_pt!=0)
2625  {
2626  // Loop over values stored at node (if any) to find
2627  // the first non-negative eqn number
2628  unsigned n_val=solid_nod_pt->variable_position_pt()->nvalue();
2629  for (unsigned i_val=0;i_val<n_val;i_val++)
2630  {
2631  int eqn_no=solid_nod_pt->
2632  variable_position_pt()->eqn_number(i_val);
2633  if (eqn_no>=0)
2634  {
2635  first_non_negative_eqn_number_plus_one=eqn_no+1;
2636  break;
2637  }
2638  }
2639  }
2640  }
2641 
2642  // Associate node with first non negative global eqn number
2643  if (first_non_negative_eqn_number_plus_one>0)
2644  {
2645  global_node_pt[first_non_negative_eqn_number_plus_one-1]=nod_pt;
2646  }
2647 
2648 
2649  // Take into account master nodes too
2650  if (dynamic_cast<RefineableElement*>(el_pt)!=0)
2651  {
2652  int n_cont_int_values=dynamic_cast<RefineableElement*>
2653  (el_pt)->ncont_interpolated_values();
2654  for (int i_cont=-1;i_cont<n_cont_int_values;i_cont++)
2655  {
2656  if (nod_pt->is_hanging(i_cont))
2657  {
2658  HangInfo* hang_pt=nod_pt->hanging_pt(i_cont);
2659  unsigned n_master=hang_pt->nmaster();
2660  for (unsigned m=0;m<n_master;m++)
2661  {
2662  Node* master_nod_pt=hang_pt->master_node_pt(m);
2663  if (!node_done[master_nod_pt])
2664  {
2665  node_done[master_nod_pt]=true;
2666 
2667  // Loop over values stored at node (if any) to find
2668  // the first non-negative eqn number
2669  unsigned first_non_negative_eqn_number_plus_one=0;
2670  unsigned n_val=master_nod_pt->nvalue();
2671  for (unsigned i_val=0;i_val<n_val;i_val++)
2672  {
2673  int eqn_no=master_nod_pt->eqn_number(i_val);
2674  if (eqn_no>=0)
2675  {
2676  first_non_negative_eqn_number_plus_one=eqn_no+1;
2677  break;
2678  }
2679  }
2680 
2681  // If we haven't found a non-negative eqn number check
2682  // eqn numbers associated with solid data (if any)
2683  if (first_non_negative_eqn_number_plus_one==0)
2684  {
2685  // If this master is a SolidNode then add its extra
2686  // eqn numbers
2687  SolidNode* master_solid_nod_pt=dynamic_cast<SolidNode*>
2688  (master_nod_pt);
2689  if (master_solid_nod_pt!=0)
2690  {
2691  // Loop over values stored at node (if any) to find
2692  // the first non-negative eqn number
2693  unsigned n_val=master_solid_nod_pt->
2694  variable_position_pt()->nvalue();
2695  for (unsigned i_val=0;i_val<n_val;i_val++)
2696  {
2697  int eqn_no=master_solid_nod_pt->
2698  variable_position_pt()->eqn_number(i_val);
2699  if (eqn_no>=0)
2700  {
2701  first_non_negative_eqn_number_plus_one=eqn_no+1;
2702  break;
2703  }
2704  }
2705  }
2706  }
2707  // Associate node with first non negative global
2708  // eqn number
2709  if (first_non_negative_eqn_number_plus_one>0)
2710  {
2711  global_node_pt[first_non_negative_eqn_number_plus_one-1]
2712  =master_nod_pt;
2713  }
2714 
2715  } // End of not-yet-done hang node
2716  }
2717  }
2718  }
2719  }
2720  } // endif for node already done
2721  }// End of loop over nodes
2722  } //End of FiniteElement
2723 
2724  // Internal data equation numbers do not need to be added since
2725  // internal data cannot be shared between distinct elements, so
2726  // internal data on locally-stored elements can never be halo.
2727  }
2728 
2729  // Set to record duplicate nodes scheduled to be killed
2730  std::set<Node*> killed_nodes;
2731 
2732  // Now loop over the other processors from highest to lowest
2733  // (i.e. if there is a duplicate between these containers
2734  // then this will use the node on the highest numbered processor)
2735  for (int iproc=n_proc-1;iproc>=0;iproc--)
2736  {
2737  // Don't have external halo elements with yourself!
2738  if (iproc!=my_rank)
2739  {
2740  // Loop over external halo elements with iproc
2741  // to remove the duplicates
2742  unsigned n_element=mesh_pt->nexternal_halo_element(iproc);
2743  for (unsigned e_ext=0;e_ext<n_element;e_ext++)
2744  {
2745  FiniteElement* finite_ext_el_pt =
2746  dynamic_cast<FiniteElement*>(mesh_pt->
2747  external_halo_element_pt(iproc,e_ext));
2748  if(finite_ext_el_pt!=0)
2749  {
2750  // Loop over nodes
2751  unsigned n_node=finite_ext_el_pt->nnode();
2752  for (unsigned j=0;j<n_node;j++)
2753  {
2754  Node* nod_pt=finite_ext_el_pt->node_pt(j);
2755 
2756  // Loop over values stored at node (if any) to find
2757  // the first non-negative eqn number
2758  unsigned first_non_negative_eqn_number_plus_one=0;
2759  unsigned n_val=nod_pt->nvalue();
2760  for (unsigned i_val=0;i_val<n_val;i_val++)
2761  {
2762  int eqn_no=nod_pt->eqn_number(i_val);
2763  if (eqn_no>=0)
2764  {
2765  first_non_negative_eqn_number_plus_one=eqn_no+1;
2766  break;
2767  }
2768  }
2769 
2770  // If we haven't found a non-negative eqn number check
2771  // eqn numbers associated with solid data (if any)
2772  if (first_non_negative_eqn_number_plus_one==0)
2773  {
2774  // Is it a solid node?
2775  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
2776  if (solid_nod_pt!=0)
2777  {
2778  // Loop over values stored at node (if any) to find
2779  // the first non-negative eqn number
2780  unsigned n_val=solid_nod_pt->variable_position_pt()->nvalue();
2781  for (unsigned i_val=0;i_val<n_val;i_val++)
2782  {
2783  int eqn_no=solid_nod_pt->
2784  variable_position_pt()->eqn_number(i_val);
2785  if (eqn_no>=0)
2786  {
2787  first_non_negative_eqn_number_plus_one=eqn_no+1;
2788  break;
2789  }
2790  }
2791  }
2792  }
2793 
2794  // Identified which node we're dealing with via first non-negative
2795  // global eqn number (if there is none, everything is pinned
2796  // and we don't give a damn...)
2797  if (first_non_negative_eqn_number_plus_one>0)
2798  {
2799  Node* existing_node_pt=
2800  global_node_pt[first_non_negative_eqn_number_plus_one-1];
2801 
2802  // Does this node already exist?
2803  if (existing_node_pt!=0)
2804  {
2805  // Record that we're about to cull one
2806  actually_removed_some_data=true;
2807 
2808  // It's a duplicate, so store the duplicated one for
2809  // later killing...
2810  Node* duplicated_node_pt=nod_pt;
2811  if (!node_done[duplicated_node_pt])
2812  {
2813 
2814  // Remove node from all boundaries
2815  std::set<unsigned>* boundaries_pt;
2816  duplicated_node_pt->
2817  get_boundaries_pt(boundaries_pt);
2818  if (boundaries_pt!=0)
2819  {
2820  Vector<unsigned> bound;
2821  unsigned nb=(*boundaries_pt).size();
2822  bound.reserve(nb);
2823  for (std::set<unsigned>::iterator it=
2824  (*boundaries_pt).begin(); it!=
2825  (*boundaries_pt).end(); it++)
2826  {
2827  bound.push_back((*it));
2828  }
2829  for (unsigned i=0;i<nb;i++)
2830  {
2831  mesh_pt->remove_boundary_node(bound[i],
2832  duplicated_node_pt);
2833  }
2834  }
2835 
2836  // Get ready to kill it
2837  killed_nodes.insert(duplicated_node_pt);
2838  unsigned i_proc=unsigned(iproc);
2839  mesh_pt->null_external_halo_node(i_proc,
2840  duplicated_node_pt);
2841  }
2842 
2843 
2844  // Note: For now we're leaving the "dangling" (no longer
2845  // accessed masters where they are; they get cleaned
2846  // up next time we delete all the external storage
2847  // for the meshes so it's a temporary "leak" only...
2848  // At some point we should probably delete them properly too
2849 
2850 #ifdef PARANOID
2851 
2852  // Check that hang status of exiting and replacement node
2853  // matches
2854  if (dynamic_cast<RefineableElement*>(finite_ext_el_pt)!=0)
2855  {
2856  int n_cont_inter_values=dynamic_cast<RefineableElement*>
2857  (finite_ext_el_pt)->ncont_interpolated_values();
2858  for (int i_cont=-1;i_cont<n_cont_inter_values;i_cont++)
2859  {
2860  unsigned n_master_orig=0;
2861  if (finite_ext_el_pt->node_pt(j)->is_hanging(i_cont))
2862  {
2863  n_master_orig=finite_ext_el_pt->node_pt(j)->
2864  hanging_pt(i_cont)->nmaster();
2865 
2866  // Temporary leak: Resolve like this:
2867  // loop over all external halo nodes and identify the
2868  // the ones that are still reached by any of the external
2869  // elements. Kill the dangling ones.
2870  }
2871  unsigned n_master_replace=0;
2872  if (existing_node_pt->is_hanging(i_cont))
2873  {
2874  n_master_replace=existing_node_pt->
2875  hanging_pt(i_cont)->nmaster();
2876  }
2877 
2878  if (n_master_orig!=n_master_replace)
2879  {
2880  std::ostringstream error_stream;
2881  error_stream << "Number of master nodes for node to be replaced, "
2882  << n_master_orig << ", doesn't match"
2883  << "those of replacement node, "
2884  << n_master_replace << " for i_cont="
2885  << i_cont << std::endl;
2886  {
2887  error_stream << "Nodal coordinates of replacement node:";
2888  unsigned ndim=existing_node_pt->ndim();
2889  for (unsigned i=0;i<ndim;i++)
2890  {
2891  error_stream << existing_node_pt->x(i) << " ";
2892  }
2893  error_stream << "\n";
2894  error_stream << "The coordinates of its "
2895  << n_master_replace << " master nodes are: \n";
2896  for (unsigned k=0;k<n_master_replace;k++)
2897  {
2898  Node* master_nod_pt=existing_node_pt->
2899  hanging_pt(i_cont)->master_node_pt(k);
2900  unsigned ndim=master_nod_pt->ndim();
2901  for (unsigned i=0;i<ndim;i++)
2902  {
2903  error_stream << master_nod_pt->x(i) << " ";
2904  }
2905  error_stream << "\n";
2906  }
2907  }
2908 
2909  {
2910  error_stream << "Nodal coordinates of node to be replaced:";
2911  unsigned ndim=finite_ext_el_pt->node_pt(j)->ndim();
2912  for (unsigned i=0;i<ndim;i++)
2913  {
2914  error_stream << finite_ext_el_pt->node_pt(j)->x(i) << " ";
2915  }
2916  error_stream << "\n";
2917  error_stream << "The coordinates of its "
2918  << n_master_orig << " master nodes are: \n";
2919  for (unsigned k=0;k<n_master_orig;k++)
2920  {
2921  Node* master_nod_pt=finite_ext_el_pt->node_pt(j)->
2922  hanging_pt(i_cont)->master_node_pt(k);
2923  unsigned ndim=master_nod_pt->ndim();
2924  for (unsigned i=0;i<ndim;i++)
2925  {
2926  error_stream << master_nod_pt->x(i) << " ";
2927  }
2928  error_stream << "\n";
2929  }
2930  }
2931 
2932 
2933  throw OomphLibError(error_stream.str(),
2934  OOMPH_CURRENT_FUNCTION,
2935  OOMPH_EXCEPTION_LOCATION);
2936  }
2937  }
2938  }
2939 #endif
2940  // ...and point to the existing one
2941  finite_ext_el_pt->node_pt(j)=existing_node_pt;
2942  }
2943  // If it doesn't add it to the list of existing ones
2944  else
2945  {
2946  global_node_pt[first_non_negative_eqn_number_plus_one-1]=
2947  nod_pt;
2948  node_done[nod_pt]=true;
2949  }
2950  }
2951 
2952 
2953  // Do the same for any master nodes of that (possibly replaced) node
2954  if (dynamic_cast<RefineableElement*>(finite_ext_el_pt)!=0)
2955  {
2956  int n_cont_inter_values=dynamic_cast<RefineableElement*>
2957  (finite_ext_el_pt)->ncont_interpolated_values();
2958  for (int i_cont=-1;i_cont<n_cont_inter_values;i_cont++)
2959  {
2960  if (finite_ext_el_pt->node_pt(j)->is_hanging(i_cont))
2961  {
2962  HangInfo* hang_pt=finite_ext_el_pt->node_pt(j)->hanging_pt(i_cont);
2963  unsigned n_master=hang_pt->nmaster();
2964  for (unsigned m=0;m<n_master;m++)
2965  {
2966  Node* master_nod_pt=hang_pt->master_node_pt(m);
2967  unsigned n_val=master_nod_pt->nvalue();
2968  unsigned first_non_negative_eqn_number_plus_one=0;
2969  for (unsigned i_val=0;i_val<n_val;i_val++)
2970  {
2971  int eqn_no=master_nod_pt->eqn_number(i_val);
2972  if (eqn_no>=0)
2973  {
2974  first_non_negative_eqn_number_plus_one=eqn_no+1;
2975  break;
2976  }
2977  }
2978 
2979  // If we haven't found a non-negative eqn number check
2980  // eqn numbers associated with solid data (if any)
2981  if (first_non_negative_eqn_number_plus_one==0)
2982  {
2983  SolidNode* solid_master_nod_pt=dynamic_cast<SolidNode*>
2984  (master_nod_pt);
2985  if (solid_master_nod_pt!=0)
2986  {
2987  // Loop over values stored at node (if any) to find
2988  // the first non-negative eqn number
2989  unsigned n_val=solid_master_nod_pt->
2990  variable_position_pt()->nvalue();
2991  for (unsigned i_val=0;i_val<n_val;i_val++)
2992  {
2993  int eqn_no=solid_master_nod_pt->
2994  variable_position_pt()->eqn_number(i_val);
2995  if (eqn_no>=0)
2996  {
2997  first_non_negative_eqn_number_plus_one=eqn_no+1;
2998  break;
2999  }
3000  }
3001  }
3002  }
3003 
3004  // Identified which node we're dealing with via
3005  // first non-negative global eqn number (if there
3006  // is none, everything is pinned and we don't give a
3007  // damn...)
3008  if (first_non_negative_eqn_number_plus_one>0)
3009  {
3010  Node* existing_node_pt=
3011  global_node_pt[
3012  first_non_negative_eqn_number_plus_one-1];
3013 
3014  // Does this node already exist?
3015  if (existing_node_pt!=0)
3016  {
3017  // Record that we're about to cull one
3018  actually_removed_some_data=true;
3019 
3020  // It's a duplicate, so store the duplicated one for
3021  // later killing...
3022  Node* duplicated_node_pt=master_nod_pt;
3023 
3024  if (!node_done[duplicated_node_pt])
3025  {
3026  // Remove node from all boundaries
3027  std::set<unsigned>* boundaries_pt;
3028  duplicated_node_pt->
3029  get_boundaries_pt(boundaries_pt);
3030  if (boundaries_pt!=0)
3031  {
3032  for (std::set<unsigned>::iterator it=
3033  (*boundaries_pt).begin(); it!=
3034  (*boundaries_pt).end(); it++)
3035  {
3036  mesh_pt->remove_boundary_node(
3037  (*it),duplicated_node_pt);
3038  }
3039  }
3040 
3041  killed_nodes.insert(duplicated_node_pt);
3042  unsigned i_proc=unsigned(iproc);
3043  mesh_pt->null_external_halo_node(i_proc,
3044  duplicated_node_pt);
3045  }
3046 
3047  // Weight of the original node
3048  double m_weight=hang_pt->master_weight(m);
3049 
3050 
3051 #ifdef PARANOID
3052  // Sanity check: setting replacement master
3053  // node for non-hanging node? Sign of really
3054  // f***ed up code.
3055  Node* tmp_nod_pt=finite_ext_el_pt->node_pt(j);
3056  if (!tmp_nod_pt->is_hanging(i_cont))
3057  {
3058  std::ostringstream error_stream;
3059  error_stream
3060  << "About to re-set master for i_cont= "
3061  << i_cont << " for external node (with proc "
3062  << iproc << " )"
3063  << tmp_nod_pt << " at ";
3064  unsigned n=tmp_nod_pt->ndim();
3065  for (unsigned jj=0;jj<n;jj++)
3066  {
3067  error_stream << tmp_nod_pt->x(jj) << " ";
3068  }
3069  error_stream
3070  << " which is not hanging --> About to die!"
3071  << "Outputting offending element into oomph-info "
3072  << "stream. \n\n";
3073  oomph_info << "\n\n";
3074  finite_ext_el_pt->output(*(oomph_info.stream_pt()));
3075  oomph_info << "\n\n";
3076  oomph_info.stream_pt()->flush();
3077  throw OomphLibError(
3078  error_stream.str(),
3079  OOMPH_CURRENT_FUNCTION,
3080  OOMPH_EXCEPTION_LOCATION);
3081  }
3082 #endif
3083 
3084 
3085  // And re-set master
3086  finite_ext_el_pt->node_pt(j)->
3087  hanging_pt(i_cont)->
3088  set_master_node_pt(m,existing_node_pt,
3089  m_weight);
3090  }
3091  // If it doesn't, add it to the list of existing ones
3092  else
3093  {
3094  global_node_pt[
3095  first_non_negative_eqn_number_plus_one-1]=
3096  master_nod_pt;
3097  node_done[master_nod_pt]=true;
3098  }
3099  }
3100  } // End of loop over master nodes
3101  } // end of hanging
3102  } // end of loop over continously interpolated variables
3103  } // end refineable element (with potentially hanging node
3104 
3105  } // end loop over nodes on external halo elements
3106 
3107  } //End of check for finite element
3108 
3109  } // end loop over external halo elements
3110  }
3111  } // end loop over processors
3112 
3113 
3114  // Now kill all the deleted nodes
3115  for (std::set<Node*>::iterator it=killed_nodes.begin();
3116  it!=killed_nodes.end();it++)
3117  {
3118  delete (*it);
3119  }
3120 
3121 
3122 // oomph_info << "Number of nonzero entries in global_node_pt: "
3123 // << global_node_pt.size() << std::endl;
3124 
3125 // // Time it...
3126 // // Taken out again by MH -- clutters up output
3127 // if (Global_timings::Doc_comprehensive_timings)
3128 // {
3129 // double t_end = TimingHelpers::timer();
3130 // oomph_info
3131 // << "Total time for Problem::remove_duplicate_data: "
3132 // << t_end-t_start << std::endl;
3133 // }
3134 
3135 }
3136 
3137 
3138 //========================================================================
3139 /// Consolidate external halo node storage by removing nulled out
3140 /// pointers in external halo and haloed schemes for all meshes.
3141 //========================================================================
3143  {
3144 
3145  //Do we have submeshes?
3146  unsigned n_mesh_loop=1;
3147  unsigned nmesh=nsub_mesh();
3148  if (nmesh>0)
3149  {
3150  n_mesh_loop=nmesh;
3151  }
3152 
3153  //Storage for number of processors and current processor
3154  int n_proc=this->communicator_pt()->nproc();
3155  int my_rank=this->communicator_pt()->my_rank();
3156 
3157  //If only one processor then return
3158  if(n_proc==1) {return;}
3159 
3160  //Loop over all (other) processors and store index of any nulled-out
3161  //external halo nodes in storage scheme.
3162 
3163  //Data to be sent to each processor
3164  Vector<int> send_n(n_proc,0);
3165 
3166  //Storage for all values to be sent to all processors
3167  Vector<int> send_data;
3168 
3169  //Start location within send_data for data to be sent to each processor
3170  Vector<int> send_displacement(n_proc,0);
3171 
3172  //Check missing ones
3173  for (int domain=0;domain<n_proc;domain++)
3174  {
3175  //Set the offset for the current processor
3176  send_displacement[domain] = send_data.size();
3177 
3178  //Don't bother to do anything if the processor in the loop is the
3179  //current processor
3180  if(domain!=my_rank)
3181  {
3182 
3183  //Deal with sub-meshes one-by-one if required
3184  Mesh* my_mesh_pt=0;
3185 
3186  //Loop over submeshes
3187  for (unsigned imesh=0;imesh<n_mesh_loop;imesh++)
3188  {
3189  if (nmesh==0)
3190  {
3191  my_mesh_pt=mesh_pt();
3192  }
3193  else
3194  {
3195  my_mesh_pt=mesh_pt(imesh);
3196  }
3197 
3198  //Make backup of external halo node pointers with this domain
3199  Vector<Node*> backup_pt(my_mesh_pt->external_halo_node_pt(domain));
3200 
3201  //How many do we have currently?
3202  unsigned nnod=backup_pt.size();
3203 
3204  //Prepare storage for updated halo nodes
3205  Vector<Node*> new_external_halo_node_pt;
3206  new_external_halo_node_pt.reserve(nnod);
3207 
3208  //Loop over external halo nodes with this domain
3209  for (unsigned j=0;j<nnod;j++)
3210  {
3211  //Get pointer to node
3212  Node* nod_pt=backup_pt[j];
3213 
3214  //Has it been nulled out?
3215  if (nod_pt==0)
3216  {
3217  //Save index of nulled out one
3218  send_data.push_back(j);
3219  }
3220  else
3221  {
3222  //Still alive: Copy across
3223  new_external_halo_node_pt.push_back(nod_pt);
3224  }
3225  }
3226 
3227  //Set new external halo node vector
3228  my_mesh_pt->set_external_halo_node_pt(domain,new_external_halo_node_pt);
3229 
3230  //End of data for this mesh
3231  send_data.push_back(-1);
3232 
3233  } // end of loop over meshes
3234 
3235  } // end skip own domain
3236 
3237  // Find the number of data added to the vector
3238  send_n[domain] = send_data.size() - send_displacement[domain];
3239 
3240  }
3241 
3242  //Storage for the number of data to be received from each processor
3243  Vector<int> receive_n(n_proc,0);
3244 
3245  //Now send numbers of data to be sent between all processors
3246  MPI_Alltoall(&send_n[0],1,MPI_INT,&receive_n[0],1,MPI_INT,
3247  this->communicator_pt()->mpi_comm());
3248 
3249 
3250  //We now prepare the data to be received
3251  //by working out the displacements from the received data
3252  Vector<int> receive_displacement(n_proc,0);
3253  int receive_data_count=0;
3254  for(int rank=0;rank<n_proc;++rank)
3255  {
3256  //Displacement is number of data received so far
3257  receive_displacement[rank] = receive_data_count;
3258  receive_data_count += receive_n[rank];
3259  }
3260 
3261  //Now resize the receive buffer for all data from all processors
3262  //Make sure that it has a size of at least one
3263  if(receive_data_count==0) {++receive_data_count;}
3264  Vector<int> receive_data(receive_data_count);
3265 
3266  //Make sure that the send buffer has size at least one
3267  //so that we don't get a segmentation fault
3268  if(send_data.size()==0) {send_data.resize(1);}
3269 
3270  //Now send the data between all the processors
3271  MPI_Alltoallv(&send_data[0],&send_n[0],&send_displacement[0],
3272  MPI_INT,
3273  &receive_data[0],&receive_n[0],
3274  &receive_displacement[0],
3275  MPI_INT,
3276  this->communicator_pt()->mpi_comm());
3277 
3278  //Now use the received data
3279  for (int send_rank=0;send_rank<n_proc;send_rank++)
3280  {
3281  //Don't bother to do anything for the processor corresponding to the
3282  //current processor or if no data were received from this processor
3283  if((send_rank != my_rank) && (receive_n[send_rank] != 0))
3284  {
3285  //Counter for the data within the large array
3286  unsigned count=receive_displacement[send_rank];
3287 
3288  //Deal with sub-meshes one-by-one if required
3289  Mesh* my_mesh_pt=0;
3290 
3291  // Loop over submeshes
3292  for (unsigned imesh=0;imesh<n_mesh_loop;imesh++)
3293  {
3294  if (nmesh==0)
3295  {
3296  my_mesh_pt=mesh_pt();
3297  }
3298  else
3299  {
3300  my_mesh_pt=mesh_pt(imesh);
3301  }
3302 
3303  //Make backup of external haloed node pointers with this domain
3304  Vector<Node*> backup_pt=my_mesh_pt->external_haloed_node_pt(send_rank);
3305 
3306  //Unpack until we reach "end of data" indicator (-1) for this mesh
3307  while (true)
3308  {
3309  //Read next entry
3310  int next_one=receive_data[count++];
3311 
3312  if (next_one==-1)
3313  {
3314  break;
3315  }
3316  else
3317  {
3318  //Null out the entry
3319  backup_pt[next_one]=0;
3320  }
3321  }
3322 
3323  //How many do we have currently?
3324  unsigned nnod=backup_pt.size();
3325 
3326  //Prepare storage for updated haloed nodes
3327  Vector<Node*> new_external_haloed_node_pt;
3328  new_external_haloed_node_pt.reserve(nnod);
3329 
3330  //Loop over external haloed nodes with this domain
3331  for (unsigned j=0;j<nnod;j++)
3332  {
3333  //Get pointer to node
3334  Node* nod_pt=backup_pt[j];
3335 
3336  //Has it been nulled out?
3337  if (nod_pt!=0)
3338  {
3339  //Still alive: Copy across
3340  new_external_haloed_node_pt.push_back(nod_pt);
3341  }
3342  }
3343 
3344  //Set new external haloed node vector
3345  my_mesh_pt->set_external_haloed_node_pt(send_rank,
3346  new_external_haloed_node_pt);
3347 
3348  }
3349  }
3350 
3351  } //End of data is received
3352  }
3353 
3354 #endif
3355 
3356 
3357 
3358 //=======================================================================
3359 /// Function that sets the values of the dofs in the object
3360 //======================================================================
3362  {
3363  const unsigned long n_dof = this->ndof();
3364 #ifdef PARANOID
3365  if(n_dof != dofs.nrow())
3366  {
3367  std::ostringstream error_stream;
3368  error_stream << "Number of degrees of freedom in vector argument "
3369  << dofs.nrow() << "\n"
3370  << "does not equal number of degrees of freedom in problem "
3371  << n_dof;
3372  throw OomphLibError(error_stream.str(),
3373  OOMPH_CURRENT_FUNCTION,
3374  OOMPH_EXCEPTION_LOCATION);
3375  }
3376 #endif
3377  for(unsigned long l=0;l<n_dof;l++)
3378  {
3379  *Dof_pt[l] = dofs[l];
3380  }
3381  }
3382 
3383  /// Set history values of dofs
3384  void Problem::set_dofs(const unsigned& t, DoubleVector& dofs)
3385  {
3386 #ifdef PARANOID
3387  if(distributed())
3388  {
3389  throw OomphLibError("Not designed for distributed problems",
3390  OOMPH_EXCEPTION_LOCATION,
3391  OOMPH_CURRENT_FUNCTION);
3392  // might work if the dofs vector is distributed in the right way...
3393  }
3394 #endif
3395 
3396  // First deal with global data
3397  unsigned Nglobal_data = nglobal_data();
3398  for(unsigned i=0; i<Nglobal_data; i++)
3399  {
3400  for(unsigned j=0, nj=Global_data_pt[i]->nvalue(); j<nj; j++)
3401  {
3402  // For each data get the equation number and copy out the value.
3403  int eqn_number = Global_data_pt[i]->eqn_number(j);
3404  if(eqn_number >= 0)
3405  {
3406  Global_data_pt[i]->set_value(t, j, dofs[eqn_number]);
3407  }
3408  }
3409  }
3410 
3411  // Next element internal data
3412  for(unsigned i=0, ni=mesh_pt()->nelement(); i<ni; i++)
3413  {
3414 
3415  GeneralisedElement* ele_pt = mesh_pt()->element_pt(i);
3416  for(unsigned j=0, nj=ele_pt->ninternal_data(); j<nj; j++)
3417  {
3418 
3419  Data* d_pt = ele_pt->internal_data_pt(j);
3420  for(unsigned k=0, nk=d_pt->nvalue(); k<nk; k++)
3421  {
3422 
3423  int eqn_number = d_pt->eqn_number(k);
3424  if(eqn_number >= 0)
3425  {
3426  d_pt->set_value(t, k, dofs[eqn_number]);
3427  }
3428  }
3429  }
3430  }
3431 
3432  // Now the nodes
3433  for(unsigned i=0, ni=mesh_pt()->nnode(); i<ni; i++)
3434  {
3435  Node* node_pt = mesh_pt()->node_pt(i);
3436  for(unsigned j=0, nj=node_pt->nvalue(); j<nj; j++)
3437  {
3438  // For each node get the equation number and copy out the value.
3439  int eqn_number = node_pt->eqn_number(j);
3440  if(eqn_number >= 0)
3441  {
3442  node_pt->set_value(t, j, dofs[eqn_number]);
3443  }
3444  }
3445  }
3446  }
3447 
3448 
3449  /// Set history values of dofs from the type of vector stored in
3450  /// problem::Dof_pt.
3451  void Problem::set_dofs(const unsigned& t, Vector<double*>& dof_pt)
3452  {
3453 #ifdef PARANOID
3454  if(distributed())
3455  {
3456  throw OomphLibError("Not implemented for distributed problems!",
3457  OOMPH_EXCEPTION_LOCATION,
3458  OOMPH_CURRENT_FUNCTION);
3459  }
3460 #endif
3461 
3462  // If we have any spine meshes I think there might be more degrees
3463  // of freedom there. I don't use them though so I'll let someone who
3464  // knows what they are doing handle it. --David Shepherd
3465 
3466  // First deal with global data
3467  unsigned Nglobal_data = nglobal_data();
3468  for(unsigned i=0; i<Nglobal_data; i++)
3469  {
3470  for(unsigned j=0, nj=Global_data_pt[i]->nvalue(); j<nj; j++)
3471  {
3472  // For each data get the equation number and copy in the value.
3473  int eqn_number = Global_data_pt[i]->eqn_number(j);
3474  if(eqn_number >= 0)
3475  {
3476  Global_data_pt[i]->set_value(t, j, *(dof_pt[eqn_number]));
3477  }
3478  }
3479  }
3480 
3481  // Now the mesh data
3482  // nodes
3483  for(unsigned i=0, ni=mesh_pt()->nnode(); i<ni; i++)
3484  {
3485  Node* node_pt = mesh_pt()->node_pt(i);
3486  for(unsigned j=0, nj=node_pt->nvalue(); j<nj; j++)
3487  {
3488  // For each node get the equation number and copy in the value.
3489  int eqn_number = node_pt->eqn_number(j);
3490  if(eqn_number >= 0)
3491  {
3492  node_pt->set_value(t, j, *(dof_pt[eqn_number]));
3493  }
3494  }
3495  }
3496 
3497  // and non-nodal data inside elements
3498  for(unsigned i=0, ni=mesh_pt()->nelement(); i<ni; i++)
3499  {
3500  GeneralisedElement* ele_pt = mesh_pt()->element_pt(i);
3501  for(unsigned j=0, nj=ele_pt->ninternal_data(); j<nj; j++)
3502  {
3503  Data* data_pt = ele_pt->internal_data_pt(j);
3504  // For each node get the equation number and copy in the value.
3505  int eqn_number = data_pt->eqn_number(j);
3506  if(eqn_number >= 0)
3507  {
3508  data_pt->set_value(t, j, *(dof_pt[eqn_number]));
3509  }
3510  }
3511  }
3512 
3513  }
3514 
3515 
3516 //===================================================================
3517 ///Function that adds the values to the dofs
3518 //==================================================================
3519  void Problem::add_to_dofs(const double &lambda,
3520  const DoubleVector &increment_dofs)
3521  {
3522  const unsigned long n_dof = this->ndof();
3523  for(unsigned long l=0;l<n_dof;l++)
3524  {
3525  *Dof_pt[l] += lambda*increment_dofs[l];
3526  }
3527  }
3528 
3529 
3530 
3531 //=========================================================================
3532 ///Return the residual vector multiplied by the inverse mass matrix
3533 ///Virtual so that it can be overloaded for mpi problems
3534 //=========================================================================
3536  {
3537  //This function does not make sense for assembly handlers other than the
3538  //default, so complain if we try to call it with another handler
3539 
3540 #ifdef PARANOID
3541  //If we are not the default, then complain
3543  {
3544  std::ostringstream error_stream;
3545  error_stream
3546  <<
3547  "The function get_inverse_mass_matrix_times_residuals() can only be\n"
3548  <<
3549  "used with the default assembly handler\n\n";
3550  throw OomphLibError(error_stream.str(),
3551  OOMPH_CURRENT_FUNCTION,
3552  OOMPH_EXCEPTION_LOCATION);
3553  }
3554 #endif
3555 
3556  //Find the number of degrees of freedom in the problem
3557  const unsigned n_dof = this->ndof();
3558 
3559  //Resize the vector
3560  LinearAlgebraDistribution dist(this->communicator_pt(),n_dof,false);
3561  Mres.build(&dist,0.0);
3562 
3563  //If we have discontinuous formulation
3564  //We can invert the mass matrix element by element
3566  {
3567  //Loop over the elements and get their residuals
3568  const unsigned n_element = Problem::mesh_pt()->nelement();
3569  Vector<double> element_Mres;
3570  for(unsigned e=0;e<n_element;e++)
3571  {
3572  //Cache the element
3573  DGElement* const elem_pt =
3574  dynamic_cast<DGElement*>(Problem::mesh_pt()->element_pt(e));
3575 
3576  //Find the elemental inverse mass matrix times residuals
3577  const unsigned n_el_dofs = elem_pt->ndof();
3578  elem_pt->get_inverse_mass_matrix_times_residuals(element_Mres);
3579 
3580  //Add contribution to global matrix
3581  for(unsigned i=0;i<n_el_dofs;i++)
3582  {
3583  Mres[elem_pt->eqn_number(i)] = element_Mres[i];
3584  }
3585  }
3586  }
3587  //Otherwise it's continous and we must invert the full
3588  //mass matrix via a global linear solve.
3589  else
3590  {
3591  //Now do the linear solve -- recycling Mass matrix if requested
3592  //If we already have the factorised mass matrix, then resolve
3594  {
3596  {
3597  oomph_info << "Not recomputing Mass Matrix " << std::endl;
3598  }
3599 
3600  //Get the residuals
3601  DoubleVector residuals(&dist,0.0);
3602  this->get_residuals(residuals);
3603 
3604  // Resolve the linear system
3606  }
3607  //Otherwise solve for the first time
3608  else
3609  {
3610  //If we wish to reuse the mass matrix, then enable resolve
3612  {
3614  {
3615  oomph_info << "Enabling resolve in explicit timestep" << std::endl;
3616  }
3618  }
3619 
3620  //Use a custom assembly handler to assemble and invert the mass matrix
3621 
3622  //Store the old assembly handler
3623  AssemblyHandler* old_assembly_handler_pt = this->assembly_handler_pt();
3624  //Set the assembly handler to the explicit timestep handler
3626 
3627  //Solve the linear system
3629  //The mass matrix has now been computed
3631 
3632  //Delete the Explicit Timestep handler
3633  delete this->assembly_handler_pt();
3634  //Reset the assembly handler to the original handler
3635  this->assembly_handler_pt() = old_assembly_handler_pt;
3636  }
3637  }
3638  }
3639 
3641  {
3642  // Loop over timesteppers: make them (temporarily) steady and store their
3643  // is_steady status.
3644  unsigned n_time_steppers = this->ntime_stepper();
3645  std::vector<bool> was_steady(n_time_steppers);
3646  for(unsigned i=0;i<n_time_steppers;i++)
3647  {
3648  was_steady[i]=time_stepper_pt(i)->is_steady();
3650  }
3651 
3652  // Calculate f using the residual/jacobian machinary.
3654 
3655  // Reset the is_steady status of all timesteppers that weren't already
3656  // steady when we came in here and reset their weights
3657  for(unsigned i=0;i<n_time_steppers;i++)
3658  {
3659  if (!was_steady[i])
3660  {
3662  }
3663  }
3664  }
3665 
3666 
3667 
3668 //================================================================
3669 /// Get the total residuals Vector for the problem
3670 //================================================================
3672  {
3673  // Three different cases; if MPI_Helpers::MPI_has_been_initialised=true
3674  // this means MPI_Helpers::init() has been called. This could happen on a
3675  // code compiled with MPI but run serially; in this instance the
3676  // get_residuals function still works on one processor.
3677  //
3678  // Secondly, if a code has been compiled with MPI, but MPI_Helpers::init()
3679  // has not been called, then MPI_Helpers::MPI_has_been_initialised=false
3680  // and the code calls...
3681  //
3682  // Thirdly, the serial version (compiled by all, but only run when compiled
3683  // with MPI if MPI_Helpers::MPI_has_been_initialised=false
3684 
3685  //Check that the residuals has the correct number of rows if it has been
3686  //setup
3687 #ifdef PARANOID
3688  if(residuals.built())
3689  {
3690  if(residuals.distribution_pt()->nrow() != this->ndof())
3691  {
3692  std::ostringstream error_stream;
3693  error_stream <<
3694  "The distribution of the residuals vector does not have the correct\n"
3695  <<
3696  "number of global rows\n";
3697 
3698  throw OomphLibError(error_stream.str(),
3699  OOMPH_CURRENT_FUNCTION,
3700  OOMPH_EXCEPTION_LOCATION);
3701  }
3702  }
3703 #endif
3704 
3705  //Find the number of rows
3706  const unsigned nrow = this->ndof();
3707 
3708  // Determine the distribution for the residuals vector
3709  // IF the vector has distribution setup then use that
3710  // ELSE determine the distribution based on the
3711  // distributed_matrix_distribution enum
3712  LinearAlgebraDistribution* dist_pt=0;
3713  if (residuals.built())
3714  {
3715  dist_pt = new LinearAlgebraDistribution(residuals.distribution_pt());
3716  }
3717  else
3718  {
3719 #ifdef OOMPH_HAS_MPI
3720  // if problem is only one one processor assemble non-distributed
3721  // distribution
3722  if (Communicator_pt->nproc() == 1)
3723  {
3724  dist_pt = new LinearAlgebraDistribution(Communicator_pt,nrow,false);
3725  }
3726  // if the problem is not distributed then assemble the jacobian with
3727  // a uniform distributed distribution
3728  else if (!Problem_has_been_distributed)
3729  {
3730  dist_pt = new LinearAlgebraDistribution(Communicator_pt,nrow,true);
3731  }
3732  // otherwise the problem is a distributed problem
3733  else
3734  {
3736  {
3738  dist_pt = new LinearAlgebraDistribution(Communicator_pt,nrow,true);
3739  break;
3742  break;
3744  LinearAlgebraDistribution* uniform_dist_pt =
3746  bool use_problem_dist = true;
3747  unsigned nproc = Communicator_pt->nproc();
3748  for (unsigned p = 0; p < nproc; p++)
3749  {
3750  if ((double)Dof_distribution_pt->nrow_local(p) >
3751  ((double)uniform_dist_pt->nrow_local(p))*1.1)
3752  {
3753  use_problem_dist = false;
3754  }
3755  }
3756  if (use_problem_dist)
3757  {
3759  }
3760  else
3761  {
3762  dist_pt = new LinearAlgebraDistribution(uniform_dist_pt);
3763  }
3764  delete uniform_dist_pt;
3765  break;
3766  }
3767  }
3768 #else
3769  dist_pt = new LinearAlgebraDistribution(Communicator_pt,nrow,false);
3770 #endif
3771  }
3772 
3773  //Locally cache pointer to assembly handler
3775 
3776  //Build and zero the residuals
3777  residuals.build(dist_pt,0.0);
3778 
3779  //Serial (or one processor case)
3780 #ifdef OOMPH_HAS_MPI
3781  if(this->communicator_pt()->nproc() == 1)
3782  {
3783 #endif // OOMPH_HAS_MPI
3784  //Loop over all the elements
3785  unsigned long Element_pt_range = Mesh_pt->nelement();
3786  for(unsigned long e=0;e<Element_pt_range;e++)
3787  {
3788  //Get the pointer to the element
3789  GeneralisedElement* elem_pt = Mesh_pt->element_pt(e);
3790  //Find number of dofs in the element
3791  unsigned n_element_dofs = assembly_handler_pt->ndof(elem_pt);
3792  //Set up an array
3793  Vector<double> element_residuals(n_element_dofs);
3794  //Fill the array
3795  assembly_handler_pt->get_residuals(elem_pt,element_residuals);
3796  //Now loop over the dofs and assign values to global Vector
3797  for(unsigned l=0;l<n_element_dofs;l++)
3798  {
3799  residuals[assembly_handler_pt->eqn_number(elem_pt,l)]
3800  += element_residuals[l];
3801  }
3802  }
3803  //Otherwise parallel case
3804 #ifdef OOMPH_HAS_MPI
3805  }
3806 else
3807  {
3808  //Store the current assembly handler
3809  AssemblyHandler* const old_assembly_handler_pt = Assembly_handler_pt;
3810  //Create a new assembly handler that only assembles the residuals
3811  Assembly_handler_pt = new ParallelResidualsHandler(old_assembly_handler_pt);
3812 
3813  //Setup memory for parallel sparse assemble
3814  //No matrix so all size zero
3815  Vector<int*> column_index;
3816  Vector<int* > row_start;
3817  Vector<double* > value;
3818  Vector<unsigned> nnz;
3819  //One set of residuals of sizer one
3820  Vector<double*> res(1);
3821 
3822  //Call the parallel sparse assemble, that should only assemble residuals
3823  parallel_sparse_assemble(dist_pt,
3824  column_index,
3825  row_start,
3826  value,
3827  nnz,
3828  res);
3829  //Fill in the residuals data
3830  residuals.set_external_values(res[0],true);
3831 
3832  //Delete new assembly handler
3833  delete Assembly_handler_pt;
3834  //Reset the assembly handler to the original
3835  Assembly_handler_pt = old_assembly_handler_pt;
3836  }
3837 #endif
3838 
3839  //Delete the distribution
3840  delete dist_pt;
3841  }
3842 
3843 //=============================================================================
3844 /// Get the fully assembled residual vector and Jacobian matrix
3845 /// in dense storage. The DoubleVector residuals returned will be
3846 /// non-distributed. If on calling this method the DoubleVector residuals is
3847 /// setup then it must be non-distributed and of the correct length.
3848 /// The matrix type DenseDoubleMatrix is not distributable and therefore
3849 /// the residual vector is also assumed to be non distributable.
3850 //=============================================================================
3852  DenseDoubleMatrix& jacobian)
3853  {
3854  // get the number of degrees of freedom
3855  unsigned n_dof=ndof();
3856 
3857 #ifdef PARANOID
3858  // PARANOID checks : if the distribution of residuals is setup then it must
3859  // must not be distributed, have the right number of rows, and the same
3860  // communicator as the problem
3861  if (residuals.built())
3862  {
3863  if (residuals.distribution_pt()->distributed())
3864  {
3865  std::ostringstream error_stream;
3866  error_stream
3867  << "If the DoubleVector residuals is setup then it must not "
3868  << "be distributed.";
3869  throw OomphLibError(error_stream.str(),
3870  OOMPH_CURRENT_FUNCTION,
3871  OOMPH_EXCEPTION_LOCATION);
3872  }
3873  if (residuals.distribution_pt()->nrow() != n_dof)
3874  {
3875  std::ostringstream error_stream;
3876  error_stream
3877  << "If the DoubleVector residuals is setup then it must have"
3878  << " the correct number of rows";
3879  throw OomphLibError(error_stream.str(),
3880  OOMPH_CURRENT_FUNCTION,
3881  OOMPH_EXCEPTION_LOCATION);
3882  }
3883  if (!(*Communicator_pt == *residuals.distribution_pt()->communicator_pt()))
3884  {
3885  std::ostringstream error_stream;
3886  error_stream
3887  << "If the DoubleVector residuals is setup then it must have"
3888  << " the same communicator as the problem.";
3889  throw OomphLibError(error_stream.str(),
3890  OOMPH_CURRENT_FUNCTION,
3891  OOMPH_EXCEPTION_LOCATION);
3892  }
3893  }
3894 #endif
3895 
3896  // set the residuals distribution if it is not setup
3897  if (!residuals.built())
3898  {
3899  LinearAlgebraDistribution dist(Communicator_pt,n_dof,false);
3900  residuals.build(&dist,0.0);
3901  }
3902  // else just zero the residuals
3903  else
3904  {
3905  residuals.initialise(0.0);
3906  }
3907 
3908  // Resize the matrices -- this cannot always be done externally
3909  // because get_jacobian exists in many different versions for
3910  // different storage formats -- resizing a CC or CR matrix doesn't
3911  // make sense.
3912 
3913  // resize the jacobian
3914  jacobian.resize(n_dof,n_dof);
3915  jacobian.initialise(0.0);
3916 
3917  //Locally cache pointer to assembly handler
3919 
3920  //Loop over all the elements
3921  unsigned long n_element = Mesh_pt->nelement();
3922  for(unsigned long e=0;e<n_element;e++)
3923  {
3924  //Get the pointer to the element
3925  GeneralisedElement* elem_pt = Mesh_pt->element_pt(e);
3926  //Find number of dofs in the element
3927  unsigned n_element_dofs = assembly_handler_pt->ndof(elem_pt);
3928  //Set up an array
3929  Vector<double> element_residuals(n_element_dofs);
3930  //Set up a matrix
3931  DenseMatrix<double> element_jacobian(n_element_dofs);
3932  //Fill the array
3933  assembly_handler_pt->get_jacobian(elem_pt,
3934  element_residuals,element_jacobian);
3935  //Now loop over the dofs and assign values to global Vector
3936  for(unsigned l=0;l<n_element_dofs;l++)
3937  {
3938  unsigned long eqn_number = assembly_handler_pt->eqn_number(elem_pt,l);
3939  residuals[eqn_number] += element_residuals[l];
3940  for(unsigned l2=0;l2<n_element_dofs;l2++)
3941  {
3942  jacobian(eqn_number ,
3943  assembly_handler_pt->eqn_number(elem_pt,l2)) +=
3944  element_jacobian(l,l2);
3945  }
3946  }
3947  }
3948  }
3949 
3950 //=============================================================================
3951 /// Return the fully-assembled Jacobian and residuals for the problem,
3952 /// in the case where the Jacobian matrix is in a distributable
3953 /// row compressed storage format.
3954 /// 1. If the distribution of the jacobian and residuals is setup then, they
3955 /// will be returned with that distribution.
3956 /// Note. the jacobian and residuals must have the same distribution.
3957 /// 2. If the distribution of the jacobian and residuals are not setup then
3958 /// their distribution will computed based on:
3959 /// Distributed_problem_matrix_distribution.
3960 //=============================================================================
3962  {
3963 
3964  // Three different cases; if MPI_Helpers::MPI_has_been_initialised=true
3965  // this means MPI_Helpers::setup() has been called. This could happen on a
3966  // code compiled with MPI but run serially; in this instance the
3967  // get_residuals function still works on one processor.
3968  //
3969  // Secondly, if a code has been compiled with MPI, but MPI_Helpers::setup()
3970  // has not been called, then MPI_Helpers::MPI_has_been_initialised=false
3971  // and the code calls...
3972  //
3973  // Thirdly, the serial version (compiled by all, but only run when compiled
3974  // with MPI if MPI_Helpers::MPI_has_been_initialised=false
3975  //
3976  // The only case where an MPI code cannot run serially at present
3977  // is one where the distribute function is used (i.e. METIS is called)
3978 
3979  //Allocate storage for the matrix entries
3980  //The generalised Vector<Vector<>> structure is required
3981  //for the most general interface to sparse_assemble() which allows
3982  //the assembly of multiple matrices at once.
3983  Vector<int* > column_index(1);
3984  Vector<int* > row_start(1);
3985  Vector<double* > value(1);
3986  Vector<unsigned> nnz(1);
3987 
3988 #ifdef PARANOID
3989  // PARANOID checks that the distribution of the jacobian matches that of the
3990  // residuals (if they are setup) and that they have the right number of rows
3991  if (residuals.built() &&
3992  jacobian.distribution_built())
3993  {
3994  if (!(*residuals.distribution_pt() == *jacobian.distribution_pt()))
3995  {
3996  std::ostringstream error_stream;
3997  error_stream << "The distribution of the residuals must "
3998  << "be the same as the distribution of the jacobian.";
3999  throw OomphLibError(error_stream.str(),
4000  OOMPH_CURRENT_FUNCTION,
4001  OOMPH_EXCEPTION_LOCATION);
4002  }
4003  if (jacobian.distribution_pt()->nrow() != this->ndof())
4004  {
4005  std::ostringstream error_stream;
4006  error_stream << "The distribution of the jacobian and residuals does not"
4007  << "have the correct number of global rows.";
4008  throw OomphLibError(error_stream.str(),
4009  OOMPH_CURRENT_FUNCTION,
4010  OOMPH_EXCEPTION_LOCATION);
4011  }
4012  }
4013  else if (residuals.built() !=
4014  jacobian.distribution_built())
4015  {
4016  std::ostringstream error_stream;
4017  error_stream << "The distribution of the jacobian and residuals must "
4018  << "both be setup or both not setup";
4019  throw OomphLibError(error_stream.str(),
4020  OOMPH_CURRENT_FUNCTION,
4021  OOMPH_EXCEPTION_LOCATION);
4022  }
4023 #endif
4024 
4025 
4026  //Allocate generalised storage format for passing to sparse_assemble()
4027  Vector<double* > res(1);
4028 
4029  // number of rows
4030  unsigned nrow = this->ndof();
4031 
4032  // determine the distribution for the jacobian.
4033  // IF the jacobian has distribution setup then use that
4034  // ELSE determine the distribution based on the
4035  // distributed_matrix_distribution enum
4036  LinearAlgebraDistribution* dist_pt=0;
4037  if (jacobian.distribution_built())
4038  {
4039  dist_pt = new LinearAlgebraDistribution(jacobian.distribution_pt());
4040  }
4041  else
4042  {
4043 #ifdef OOMPH_HAS_MPI
4044  // if problem is only one one processor
4045  if (Communicator_pt->nproc() == 1)
4046  {
4047  dist_pt = new LinearAlgebraDistribution(Communicator_pt,nrow,false);
4048  }
4049  // if the problem is not distributed then assemble the jacobian with
4050  // a uniform distributed distribution
4051  else if (!Problem_has_been_distributed)
4052  {
4053  dist_pt = new LinearAlgebraDistribution(Communicator_pt,nrow,true);
4054  }
4055  // otherwise the problem is a distributed problem
4056  else
4057  {
4059  {
4061  dist_pt = new LinearAlgebraDistribution(Communicator_pt,nrow,true);
4062  break;
4065  break;
4067  LinearAlgebraDistribution* uniform_dist_pt =
4069  bool use_problem_dist = true;
4070  unsigned nproc = Communicator_pt->nproc();
4071  for (unsigned p = 0; p < nproc; p++)
4072  {
4073  if ((double)Dof_distribution_pt->nrow_local(p) >
4074  ((double)uniform_dist_pt->nrow_local(p))*1.1)
4075  {
4076  use_problem_dist = false;
4077  }
4078  }
4079  if (use_problem_dist)
4080  {
4082  }
4083  else
4084  {
4085  dist_pt = new LinearAlgebraDistribution(uniform_dist_pt);
4086  }
4087  delete uniform_dist_pt;
4088  break;
4089  }
4090  }
4091 #else
4092  dist_pt = new LinearAlgebraDistribution(Communicator_pt,nrow,false);
4093 #endif
4094  }
4095 
4096 
4097  //The matrix is in compressed row format
4098  bool compressed_row_flag=true;
4099 
4100 #ifdef OOMPH_HAS_MPI
4101  //
4102  if (Communicator_pt->nproc() == 1)
4103  {
4104 #endif
4106  row_start,
4107  value,
4108  nnz,
4109  res,
4110  compressed_row_flag);
4111  jacobian.build(dist_pt);
4112  jacobian.build_without_copy(dist_pt->nrow(),nnz[0],
4113  value[0],column_index[0],row_start[0]);
4114  residuals.build(dist_pt,0.0);
4115  residuals.set_external_values(res[0],true);
4116 #ifdef OOMPH_HAS_MPI
4117  }
4118  else
4119  {
4120  if (dist_pt->distributed())
4121  {
4122  parallel_sparse_assemble(dist_pt,
4123  column_index,
4124  row_start,
4125  value,
4126  nnz,
4127  res);
4128  jacobian.build(dist_pt);
4129  jacobian.build_without_copy(dist_pt->nrow(),nnz[0],
4130  value[0],column_index[0],
4131  row_start[0]);
4132  residuals.build(dist_pt,0.0);
4133  residuals.set_external_values(res[0],true);
4134  }
4135  else
4136  {
4137  LinearAlgebraDistribution* temp_dist_pt =
4138  new LinearAlgebraDistribution(Communicator_pt,dist_pt->nrow(),true);
4139  parallel_sparse_assemble(temp_dist_pt,
4140  column_index,
4141  row_start,
4142  value,
4143  nnz,
4144  res);
4145  jacobian.build(temp_dist_pt);
4146  jacobian.build_without_copy(dist_pt->nrow(),nnz[0],
4147  value[0],column_index[0],
4148  row_start[0]);
4149  jacobian.redistribute(dist_pt);
4150  residuals.build(temp_dist_pt,0.0);
4151  residuals.set_external_values(res[0],true);
4152  residuals.redistribute(dist_pt);
4153  delete temp_dist_pt;
4154  }
4155  }
4156 #endif
4157 
4158  // clean up dist_pt and residuals_vector pt
4159  delete dist_pt;
4160  }
4161 
4162 //=============================================================================
4163 /// Return the fully-assembled Jacobian and residuals for the problem,
4164 /// in the case when the jacobian matrix is in column-compressed storage
4165 /// format.
4166 //=============================================================================
4168  {
4169  // Three different cases; if MPI_Helpers::MPI_has_been_initialised=true
4170  // this means MPI_Helpers::setup() has been called. This could happen on a
4171  // code compiled with MPI but run serially; in this instance the
4172  // get_residuals function still works on one processor.
4173  //
4174  // Secondly, if a code has been compiled with MPI, but MPI_Helpers::setup()
4175  // has not been called, then MPI_Helpers::MPI_has_been_initialised=false
4176  // and the code calls...
4177  //
4178  // Thirdly, the serial version (compiled by all, but only run when compiled
4179  // with MPI if MPI_Helpers::MPI_has_been_5Binitialised=false
4180  //
4181  // The only case where an MPI code cannot run serially at present
4182  // is one where the distribute function is used (i.e. METIS is called)
4183 
4184  // get the number of degrees of freedom
4185  unsigned n_dof=ndof();
4186 
4187 #ifdef PARANOID
4188  // PARANOID checks : if the distribution of residuals is setup then it must
4189  // must not be distributed, have the right number of rows, and the same
4190  // communicator as the problem
4191  if (residuals.built())
4192  {
4193  if (residuals.distribution_pt()->distributed())
4194  {
4195  std::ostringstream error_stream;
4196  error_stream
4197  << "If the DoubleVector residuals is setup then it must not "
4198  << "be distributed.";
4199  throw OomphLibError(error_stream.str(),
4200  OOMPH_CURRENT_FUNCTION,
4201  OOMPH_EXCEPTION_LOCATION);
4202  }
4203  if (residuals.distribution_pt()->nrow() != n_dof)
4204  {
4205  std::ostringstream error_stream;
4206  error_stream
4207  << "If the DoubleVector residuals is setup then it must have"
4208  << " the correct number of rows";
4209  throw OomphLibError(error_stream.str(),
4210  OOMPH_CURRENT_FUNCTION,
4211  OOMPH_EXCEPTION_LOCATION);
4212  }
4213  if (!(*Communicator_pt == *residuals.distribution_pt()->communicator_pt()))
4214  {
4215  std::ostringstream error_stream;
4216  error_stream
4217  << "If the DoubleVector residuals is setup then it must have"
4218  << " the same communicator as the problem.";
4219  throw OomphLibError(error_stream.str(),
4220  OOMPH_CURRENT_FUNCTION,
4221  OOMPH_EXCEPTION_LOCATION);
4222  }
4223  }
4224 #endif
4225 
4226  //Allocate storage for the matrix entries
4227  //The generalised Vector<Vector<>> structure is required
4228  //for the most general interface to sparse_assemble() which allows
4229  //the assembly of multiple matrices at once.
4230  Vector<int* > row_index(1);
4231  Vector<int* > column_start(1);
4232  Vector<double* > value(1);
4233 
4234  //Allocate generalised storage format for passing to sparse_assemble()
4235  Vector<double* > res(1);
4236 
4237  // allocate storage for the number of non-zeros in each matrix
4238  Vector<unsigned> nnz(1);
4239 
4240  //The matrix is in compressed column format
4241  bool compressed_row_flag=false;
4242 
4243  // get the distribution for the residuals
4244  LinearAlgebraDistribution* dist_pt;
4245  if (!residuals.built())
4246  {
4247  dist_pt
4248  = new LinearAlgebraDistribution(Communicator_pt,this->ndof(),false);
4249  }
4250  else
4251  {
4252  dist_pt = new LinearAlgebraDistribution(residuals.distribution_pt());
4253  }
4254 
4255 #ifdef OOMPH_HAS_MPI
4256  if (communicator_pt()->nproc() == 1)
4257  {
4258 #endif
4260  column_start,
4261  value,
4262  nnz,
4263  res,
4264  compressed_row_flag);
4265  jacobian.build_without_copy(value[0],row_index[0],column_start[0],nnz[0],
4266  n_dof,n_dof);
4267  residuals.build(dist_pt,0.0);
4268  residuals.set_external_values(res[0],true);
4269 #ifdef OOMPH_HAS_MPI
4270  }
4271  else
4272  {
4273  std::ostringstream error_stream;
4274  error_stream
4275  << "Cannot assemble a CCDoubleMatrix Jacobian on more "
4276  << "than one processor.";
4277  throw OomphLibError(error_stream.str(),
4278  OOMPH_CURRENT_FUNCTION,
4279  OOMPH_EXCEPTION_LOCATION);
4280  }
4281 #endif
4282 
4283  // clean up
4284  delete dist_pt;
4285  }
4286 
4287 
4288 //===================================================================
4289 /// \short Set all pinned values to zero.
4290 /// Used to set boundary conditions to be homogeneous in the copy
4291 /// of the problem used in adaptive bifurcation tracking
4292 /// (ALH: TEMPORARY HACK, WILL BE FIXED)
4293 //==================================================================
4295 {
4296  //NOTE THIS DOES NOT ZERO ANY SPINE DATA, but otherwise everything else
4297  //should be zeroed
4298 
4299  //Zero any pinned global Data
4300  const unsigned n_global_data = nglobal_data();
4301  for(unsigned i=0;i<n_global_data;i++)
4302  {
4303  Data* const local_data_pt = Global_data_pt[i];
4304  const unsigned n_value = local_data_pt->nvalue();
4305  for(unsigned j=0;j<n_value;j++)
4306  {
4307  //If the data value is pinned set the value to zero
4308  if(local_data_pt->is_pinned(j)) {local_data_pt->set_value(j,0.0);}
4309  }
4310  }
4311 
4312  // Loop over the submeshes:
4313  const unsigned n_sub_mesh=Sub_mesh_pt.size();
4314  if(n_sub_mesh==0)
4315  {
4316  //Loop over the nodes in the element
4317  const unsigned n_node = Mesh_pt->nnode();
4318  for(unsigned n=0;n<n_node;n++)
4319  {
4320  Node* const local_node_pt = Mesh_pt->node_pt(n);
4321  const unsigned n_value = local_node_pt->nvalue();
4322  for(unsigned j=0;j<n_value;j++)
4323  {
4324  //If the data value is pinned set the value to zero
4325  if(local_node_pt->is_pinned(j)) {local_node_pt->set_value(j,0.0);}
4326  }
4327 
4328  //Try to cast to a solid node
4329  SolidNode* const local_solid_node_pt =
4330  dynamic_cast<SolidNode*>(local_node_pt);
4331  //If we are successful
4332  if(local_solid_node_pt)
4333  {
4334  //Find the dimension of the node
4335  const unsigned n_dim = local_solid_node_pt->ndim();
4336  //Find number of positions
4337  const unsigned n_position_type = local_solid_node_pt->nposition_type();
4338 
4339  for(unsigned k=0;k<n_position_type;k++)
4340  {
4341  for(unsigned i=0;i<n_dim;i++)
4342  {
4343  //If the generalised position is pinned,
4344  //set the value to zero
4345  if(local_solid_node_pt->position_is_pinned(k,i))
4346  {
4347  local_solid_node_pt->x_gen(k,i) = 0.0;
4348  }
4349  }
4350  }
4351  }
4352  }
4353 
4354  //Now loop over the element's and zero the internal data
4355  const unsigned n_element = Mesh_pt->nelement();
4356  for(unsigned e=0;e<n_element;e++)
4357  {
4358  GeneralisedElement* const local_element_pt = Mesh_pt->element_pt(e);
4359  const unsigned n_internal = local_element_pt->ninternal_data();
4360  for(unsigned i=0;i<n_internal;i++)
4361  {
4362  Data* const local_data_pt = local_element_pt->internal_data_pt(i);
4363  const unsigned n_value = local_data_pt->nvalue();
4364  for(unsigned j=0;j<n_value;j++)
4365  {
4366  //If the data value is pinned set the value to zero
4367  if(local_data_pt->is_pinned(j)) {local_data_pt->set_value(j,0.0);}
4368  }
4369  }
4370  } //End of loop over elements
4371  }
4372  else
4373  {
4374  //Alternatively loop over all sub meshes
4375  for (unsigned m=0;m<n_sub_mesh;m++)
4376  {
4377  //Loop over the nodes in the element
4378  const unsigned n_node = Sub_mesh_pt[m]->nnode();
4379  for(unsigned n=0;n<n_node;n++)
4380  {
4381  Node* const local_node_pt = Sub_mesh_pt[m]->node_pt(n);
4382  const unsigned n_value = local_node_pt->nvalue();
4383  for(unsigned j=0;j<n_value;j++)
4384  {
4385  //If the data value is pinned set the value to zero
4386  if(local_node_pt->is_pinned(j)) {local_node_pt->set_value(j,0.0);}
4387  }
4388 
4389  //Try to cast to a solid node
4390  SolidNode* const local_solid_node_pt =
4391  dynamic_cast<SolidNode*>(local_node_pt);
4392  //If we are successful
4393  if(local_solid_node_pt)
4394  {
4395  //Find the dimension of the node
4396  const unsigned n_dim = local_solid_node_pt->ndim();
4397  //Find number of positions
4398  const unsigned n_position_type =
4399  local_solid_node_pt->nposition_type();
4400 
4401  for(unsigned k=0;k<n_position_type;k++)
4402  {
4403  for(unsigned i=0;i<n_dim;i++)
4404  {
4405  //If the generalised position is pinned,
4406  //set the value to zero
4407  if(local_solid_node_pt->position_is_pinned(k,i))
4408  {
4409  local_solid_node_pt->x_gen(k,i) = 0.0;
4410  }
4411  }
4412  }
4413  }
4414  }
4415 
4416  //Now loop over the element's and zero the internal data
4417  const unsigned n_element = Sub_mesh_pt[m]->nelement();
4418  for(unsigned e=0;e<n_element;e++)
4419  {
4420  GeneralisedElement* const local_element_pt =
4421  Sub_mesh_pt[m]->element_pt(e);
4422  const unsigned n_internal = local_element_pt->ninternal_data();
4423  for(unsigned i=0;i<n_internal;i++)
4424  {
4425  Data* const local_data_pt = local_element_pt->internal_data_pt(i);
4426  const unsigned n_value = local_data_pt->nvalue();
4427  for(unsigned j=0;j<n_value;j++)
4428  {
4429  //If the data value is pinned set the value to zero
4430  if(local_data_pt->is_pinned(j)) {local_data_pt->set_value(j,0.0);}
4431  }
4432  }
4433  } //End of loop over elements
4434  }
4435  }
4436 }
4437 
4438 
4439 
4440 
4441 //=====================================================================
4442 /// This is a (private) helper function that is used to assemble system
4443 /// matrices in compressed row or column format
4444 /// and compute residual vectors.
4445 /// The default action is to assemble the jacobian matrix and
4446 /// residuals for the Newton method. The action can be
4447 /// overloaded at an elemental level by changing the default
4448 /// behaviour of the function Element::get_all_vectors_and_matrices().
4449 /// column_or_row_index: Column [or row] index of given entry
4450 /// row_or_column_start: Index of first entry for given row [or column]
4451 /// value : Vector of nonzero entries
4452 /// residuals : Residual vector
4453 /// compressed_row_flag: Bool flag to indicate if storage format is
4454 /// compressed row [if false interpretation of
4455 /// arguments is as stated in square brackets].
4456 /// We provide four different assembly methods, each with different
4457 /// memory requirements/execution speeds. The method is set by
4458 /// the public flag Problem::Sparse_assembly_method.
4459 //=====================================================================
4461  Vector<int* > &column_or_row_index,
4462  Vector<int* > &row_or_column_start,
4463  Vector<double* > &value,
4464  Vector<unsigned > &nnz,
4465  Vector<double* > &residuals,
4466  bool compressed_row_flag)
4467 {
4468 
4469  // Choose the actual method
4470  switch(Sparse_assembly_method)
4471  {
4472 
4474 
4476  column_or_row_index,
4477  row_or_column_start,
4478  value,
4479  nnz,
4480  residuals,
4481  compressed_row_flag);
4482 
4483  break;
4484 
4486 
4488  column_or_row_index,
4489  row_or_column_start,
4490  value,
4491  nnz,
4492  residuals,
4493  compressed_row_flag);
4494 
4495  break;
4496 
4498 
4500  column_or_row_index,
4501  row_or_column_start,
4502  value,
4503  nnz,
4504  residuals,
4505  compressed_row_flag);
4506 
4507  break;
4508 
4510 
4512  column_or_row_index,
4513  row_or_column_start,
4514  value,
4515  nnz,
4516  residuals,
4517  compressed_row_flag);
4518 
4519  break;
4520 
4522 
4524  column_or_row_index,
4525  row_or_column_start,
4526  value,
4527  nnz,
4528  residuals,
4529  compressed_row_flag);
4530 
4531  break;
4532 
4533  default:
4534 
4535  std::ostringstream error_stream;
4536  error_stream
4537  << "Error: Incorrect value for Problem::Sparse_assembly_method"
4538  << Sparse_assembly_method << std::endl
4539  << "It should be one of the enumeration Problem::Assembly_method"
4540  << std::endl;
4541  throw OomphLibError(error_stream.str(),
4542  OOMPH_CURRENT_FUNCTION,
4543  OOMPH_EXCEPTION_LOCATION);
4544  }
4545 }
4546 
4547 
4548 
4549 
4550 
4551 //=====================================================================
4552 /// This is a (private) helper function that is used to assemble system
4553 /// matrices in compressed row or column format
4554 /// and compute residual vectors, using maps
4555 /// The default action is to assemble the jacobian matrix and
4556 /// residuals for the Newton method. The action can be
4557 /// overloaded at an elemental level by chaging the default
4558 /// behaviour of the function Element::get_all_vectors_and_matrices().
4559 /// column_or_row_index: Column [or row] index of given entry
4560 /// row_or_column_start: Index of first entry for given row [or column]
4561 /// value : Vector of nonzero entries
4562 /// residuals : Residual vector
4563 /// compressed_row_flag: Bool flag to indicate if storage format is
4564 /// compressed row [if false interpretation of
4565 /// arguments is as stated in square brackets].
4566 //=====================================================================
4568  Vector<int* > &column_or_row_index,
4569  Vector<int* > &row_or_column_start,
4570  Vector<double* > &value,
4571  Vector<unsigned> &nnz,
4572  Vector<double* > &residuals,
4573  bool compressed_row_flag)
4574 {
4575  //Total number of elements
4576  const unsigned long n_elements = mesh_pt()->nelement();
4577 
4578  // Default range of elements for distributed problems
4579  unsigned long el_lo=0;
4580  unsigned long el_hi=n_elements-1;
4581 
4582 #ifdef OOMPH_HAS_MPI
4583  // Otherwise just loop over a fraction of the elements
4584  // (This will either have been initialised in
4585  // Problem::set_default_first_and_last_element_for_assembly() or
4586  // will have been re-assigned during a previous assembly loop
4587  // Note that following the re-assignment only the entries
4588  // for the current processor are relevant.
4590  {
4591  el_lo=First_el_for_assembly[Communicator_pt->my_rank()];
4592  el_hi=Last_el_plus_one_for_assembly[Communicator_pt->my_rank()]-1;
4593  }
4594 #endif
4595 
4596  // number of dofs
4597  unsigned ndof = this->ndof();
4598 
4599  //Find the number of vectors to be assembled
4600  const unsigned n_vector = residuals.size();
4601 
4602  //Find the number of matrices to be assembled
4603  const unsigned n_matrix = column_or_row_index.size();
4604 
4605  //Locally cache pointer to assembly handler
4607 
4608 #ifdef OOMPH_HAS_MPI
4609  bool doing_residuals=false;
4610  if (dynamic_cast<ParallelResidualsHandler*>(Assembly_handler_pt)!=0)
4611  {
4612  doing_residuals=true;
4613  }
4614 #endif
4615 
4616 //Error check dimensions
4617 #ifdef PARANOID
4618  if(row_or_column_start.size() != n_matrix)
4619  {
4620  std::ostringstream error_stream;
4621  error_stream
4622  << "Error: " << std::endl
4623  << "row_or_column_start.size() " << row_or_column_start.size()
4624  << " does not equal "
4625  << "column_or_row_index.size() "
4626  << column_or_row_index.size() << std::endl;
4627  throw OomphLibError(
4628  error_stream.str(),
4629  OOMPH_CURRENT_FUNCTION,
4630  OOMPH_EXCEPTION_LOCATION);
4631  }
4632 
4633  if(value.size() != n_matrix)
4634  {
4635  std::ostringstream error_stream;
4636  error_stream
4637  << "Error in Problem::sparse_assemble_row_or_column_compressed "
4638  << std::endl
4639  << "value.size() " << value.size() << " does not equal "
4640  << "column_or_row_index.size() "
4641  << column_or_row_index.size() << std::endl<< std::endl
4642  << std::endl;
4643  throw OomphLibError(
4644  error_stream.str(),
4645  OOMPH_CURRENT_FUNCTION,
4646  OOMPH_EXCEPTION_LOCATION);
4647  }
4648 #endif
4649 
4650 
4651 //The idea behind this sparse assembly routine is to use a vector of
4652 //maps for the entries in each row or column of the complete matrix.
4653 //The key for each map is the global row or column number and
4654 //the default comparison operator for integers means that each map
4655 //is ordered by the global row or column number. Thus, we need not
4656 //sort the maps, that happens at each insertion of a new entry. The
4657 //price we pay is that for large maps, inseration is not a
4658 //cheap operation. Hash maps can be used to increase the speed, but then
4659 //the ordering is lost and we would have to sort anyway. The solution if
4660 //speed is required is to use lists, see below.
4661 
4662 
4663 //Set up a vector of vectors of maps of entries of each matrix,
4664 //indexed by either the column or row. The entries of the vector for
4665 //each matrix correspond to all the rows or columns of that matrix.
4666 //The use of the map storage
4667 //scheme, with its implicit ordering on the first index, gives
4668 //a sparse ordered list of the entries in the given row or column.
4669  Vector<Vector<std::map<unsigned,double> > > matrix_data_map(n_matrix);
4670  //Loop over the number of matrices being assembled and resize
4671  //each vector of maps to the number of rows or columns of the matrix
4672  for(unsigned m=0;m<n_matrix;m++) {matrix_data_map[m].resize(ndof);}
4673 
4674  //Resize the residuals vectors
4675  for(unsigned v=0;v<n_vector;v++)
4676  {
4677  residuals[v] = new double[ndof];
4678  for (unsigned i = 0; i < ndof; i++)
4679  {
4680  residuals[v][i] = 0;
4681  }
4682  }
4683 
4684 
4685 #ifdef OOMPH_HAS_MPI
4686 
4687 
4688  // Storage for assembly time for elements
4689  double t_assemble_start=0.0;
4690 
4691  // Storage for assembly times
4692  if ((!doing_residuals)&&
4694  {
4695  Elemental_assembly_time.resize(n_elements);
4696  }
4697 
4698 #endif
4699 
4700  //----------------Assemble and populate the maps-------------------------
4701  {
4702  //Allocate local storage for the element's contribution to the
4703  //residuals vectors and system matrices of the size of the maximum
4704  //number of dofs in any element.
4705  //This means that the storage is only allocated (and deleted) once
4706  Vector<Vector<double> > el_residuals(n_vector);
4707  Vector<DenseMatrix<double> > el_jacobian(n_matrix);
4708 
4709  //Loop over the elements for this processor
4710  for(unsigned long e=el_lo;e<=el_hi;e++)
4711  {
4712 
4713 #ifdef OOMPH_HAS_MPI
4714  // Time it?
4715  if ((!doing_residuals)&&
4717  {
4718  t_assemble_start=TimingHelpers::timer();
4719  }
4720 #endif
4721 
4722  //Get the pointer to the element
4723  GeneralisedElement* elem_pt = mesh_pt()->element_pt(e);
4724 
4725 #ifdef OOMPH_HAS_MPI
4726  //Ignore halo elements
4727  if (!elem_pt->is_halo())
4728  {
4729 #endif
4730 
4731  //Find number of degrees of freedom in the element
4732  const unsigned nvar = assembly_handler_pt->ndof(elem_pt);
4733 
4734  //Resize the storage for elemental jacobian and residuals
4735  for(unsigned v=0;v<n_vector;v++) {el_residuals[v].resize(nvar);}
4736  for(unsigned m=0;m<n_matrix;m++) {el_jacobian[m].resize(nvar);}
4737 
4738  //Now get the residuals and jacobian for the element
4739  assembly_handler_pt->
4740  get_all_vectors_and_matrices(elem_pt,el_residuals, el_jacobian);
4741 
4742  //---------------Insert the values into the maps--------------
4743 
4744  //Loop over the first index of local variables
4745  for(unsigned i=0;i<nvar;i++)
4746  {
4747  //Get the local equation number
4748  unsigned eqn_number
4749  = assembly_handler_pt->eqn_number(elem_pt,i);
4750 
4751  //Add the contribution to the residuals
4752  for(unsigned v=0;v<n_vector;v++)
4753  {
4754  //Fill in each residuals vector
4755  residuals[v][eqn_number] += el_residuals[v][i];
4756  }
4757 
4758  //Now loop over the other index
4759  for(unsigned j=0;j<nvar;j++)
4760  {
4761  //Get the number of the unknown
4762  unsigned unknown = assembly_handler_pt->eqn_number(elem_pt,j);
4763 
4764  //Loop over the matrices
4765  for(unsigned m=0;m<n_matrix;m++)
4766  {
4767  //Get the value of the matrix at this point
4768  double value = el_jacobian[m](i,j);
4769  //Only bother to add to the map if it's non-zero
4770  if(std::fabs(value) > Numerical_zero_for_sparse_assembly)
4771  {
4772  //If it's compressed row storage, then our vector of maps
4773  //is indexed by row (equation number)
4774  if(compressed_row_flag)
4775  {
4776  //Add the data into the map using the unknown as the map key
4777  matrix_data_map[m][eqn_number][unknown] += value;
4778  }
4779  //Otherwise it's compressed column storage and our vector is
4780  //indexed by column (the unknown)
4781  else
4782  {
4783  //Add the data into the map using the eqn_numbe as the map key
4784  matrix_data_map[m][unknown][eqn_number] += value;
4785  }
4786  }
4787  } //End of loop over matrices
4788  }
4789  }
4790 
4791 #ifdef OOMPH_HAS_MPI
4792  } // endif halo element
4793 #endif
4794 
4795 
4796 #ifdef OOMPH_HAS_MPI
4797  // Time it?
4798  if ((!doing_residuals)&&
4800  {
4801  Elemental_assembly_time[e]=TimingHelpers::timer()-t_assemble_start;
4802  }
4803 #endif
4804 
4805  } //End of loop over the elements
4806 
4807  } //End of map assembly
4808 
4809 
4810 #ifdef OOMPH_HAS_MPI
4811 
4812  // Postprocess timing information and re-allocate distribution of
4813  // elements during subsequent assemblies.
4814  if ((!doing_residuals)&&
4817  {
4819  }
4820 
4821  // We have determined load balancing for current setup.
4822  // This can remain the same until assign_eqn_numbers() is called
4823  // again -- the flag is re-set to true there.
4824  if ((!doing_residuals)&&
4826  {
4828  }
4829 
4830 #endif
4831 
4832 
4833  //-----------Finally we need to convert the beautiful map storage scheme
4834  //------------------------to the containers required by SuperLU
4835 
4836  //Loop over the number of matrices
4837  for(unsigned m=0;m<n_matrix;m++)
4838  {
4839  //Set the number of rows or columns
4840  row_or_column_start[m] = new int[ndof+1];
4841  //Counter for the total number of entries in the storage scheme
4842  unsigned long entry_count=0;
4843  row_or_column_start[m][0] = entry_count;
4844 
4845  // first we compute the number of non-zeros
4846  nnz[m] = 0;
4847  for(unsigned long i_global=0;i_global<ndof;i_global++)
4848  {
4849  nnz[m] += matrix_data_map[m][i_global].size();
4850  }
4851 
4852  // and then resize the storage
4853  column_or_row_index[m] = new int[nnz[m]];
4854  value[m] = new double[nnz[m]];
4855 
4856  //Now we merely loop over the number of rows or columns
4857  for(unsigned long i_global=0;i_global<ndof;i_global++)
4858  {
4859  //Start index for the present row
4860  row_or_column_start[m][i_global] = entry_count;
4861  //If there are no entries in the map then skip the rest of the loop
4862  if(matrix_data_map[m][i_global].empty()) {continue;}
4863 
4864  //Loop over all the entries in the map corresponding to the given
4865  //row or column. It will be ordered
4866 
4867  for(std::map<unsigned,double>::iterator
4868  it = matrix_data_map[m][i_global].begin();
4869  it!=matrix_data_map[m][i_global].end();++it)
4870  {
4871  //The first value is the column or row index
4872  column_or_row_index[m][entry_count] = it->first;
4873  //The second value is the actual data value
4874  value[m][entry_count] = it->second;
4875  //Increase the value of the counter
4876  entry_count++;
4877  }
4878  }
4879 
4880  //Final entry in the row/column start vector
4881  row_or_column_start[m][ndof] = entry_count;
4882  } //End of the loop over the matrices
4883 
4885  {
4886  oomph_info << "Pausing at end of sparse assembly." << std::endl;
4887  pause("Check memory usage now.");
4888  }
4889 }
4890 
4891 
4892 
4893 
4894 
4895 
4896 //=====================================================================
4897 /// This is a (private) helper function that is used to assemble system
4898 /// matrices in compressed row or column format
4899 /// and compute residual vectors using lists
4900 /// The default action is to assemble the jacobian matrix and
4901 /// residuals for the Newton method. The action can be
4902 /// overloaded at an elemental level by chaging the default
4903 /// behaviour of the function Element::get_all_vectors_and_matrices().
4904 /// column_or_row_index: Column [or row] index of given entry
4905 /// row_or_column_start: Index of first entry for given row [or column]
4906 /// value : Vector of nonzero entries
4907 /// residuals : Residual vector
4908 /// compressed_row_flag: Bool flag to indicate if storage format is
4909 /// compressed row [if false interpretation of
4910 /// arguments is as stated in square brackets].
4911 //=====================================================================
4913  Vector<int* > &column_or_row_index,
4914  Vector<int* > &row_or_column_start,
4915  Vector<double* > &value,
4916  Vector<unsigned> &nnz,
4917  Vector<double* > &residuals,
4918  bool compressed_row_flag)
4919 {
4920  //Total number of elements
4921  const unsigned long n_elements = mesh_pt()->nelement();
4922 
4923  // Default range of elements for distributed problems
4924  unsigned long el_lo=0;
4925  unsigned long el_hi=n_elements-1;
4926 
4927 #ifdef OOMPH_HAS_MPI
4928  // Otherwise just loop over a fraction of the elements
4929  // (This will either have been initialised in
4930  // Problem::set_default_first_and_last_element_for_assembly() or
4931  // will have been re-assigned during a previous assembly loop
4932  // Note that following the re-assignment only the entries
4933  // for the current processor are relevant.
4935  {
4936  el_lo=First_el_for_assembly[Communicator_pt->my_rank()];
4937  el_hi=Last_el_plus_one_for_assembly[Communicator_pt->my_rank()]-1;
4938  }
4939 #endif
4940 
4941  // number of dofs
4942  unsigned ndof = this->ndof();
4943 
4944  //Find the number of vectors to be assembled
4945  const unsigned n_vector = residuals.size();
4946 
4947  //Find the number of matrices to be assembled
4948  const unsigned n_matrix = column_or_row_index.size();
4949 
4950  //Locally cache pointer to assembly handler
4952 
4953 #ifdef OOMPH_HAS_MPI
4954  bool doing_residuals=false;
4955  if (dynamic_cast<ParallelResidualsHandler*>(Assembly_handler_pt)!=0)
4956  {
4957  doing_residuals=true;
4958  }
4959 #endif
4960 
4961 //Error check dimensions
4962 #ifdef PARANOID
4963  if(row_or_column_start.size() != n_matrix)
4964  {
4965  std::ostringstream error_stream;
4966  error_stream
4967  << "Error: " << std::endl
4968  << "row_or_column_start.size() " << row_or_column_start.size()
4969  << " does not equal "
4970  << "column_or_row_index.size() "
4971  << column_or_row_index.size() << std::endl;
4972  throw OomphLibError(
4973  error_stream.str(),
4974  OOMPH_CURRENT_FUNCTION,
4975  OOMPH_EXCEPTION_LOCATION);
4976  }
4977 
4978  if(value.size() != n_matrix)
4979  {
4980  std::ostringstream error_stream;
4981  error_stream
4982  << "Error in Problem::sparse_assemble_row_or_column_compressed "
4983  << std::endl
4984  << "value.size() " << value.size() << " does not equal "
4985  << "column_or_row_index.size() "
4986  << column_or_row_index.size() << std::endl<< std::endl
4987  << std::endl;
4988  throw OomphLibError(
4989  error_stream.str(),
4990  OOMPH_CURRENT_FUNCTION,
4991  OOMPH_EXCEPTION_LOCATION);
4992  }
4993 #endif
4994 
4995 //The idea behind this sparse assembly routine is to use a vector of
4996 //lists for the entries in each row or column of the complete matrix.
4997 //The lists contain pairs of entries (global row/column number, value).
4998 //All non-zero contributions from each element are added to the lists.
4999 //We then sort each list by global row/column number and then combine
5000 //the entries corresponding to each row/column before adding to the
5001 //vectors column_or_row_index and value.
5002 
5003 //Note the trade off for "fast assembly" is that we will require
5004 //more memory during the assembly phase. Then again, if we can
5005 //only just assemble the sparse matrix, we're in real trouble.
5006 
5007 //Set up a vector of lists of paired entries of
5008 //(row/column index, jacobian matrix entry).
5009 //The entries of the vector correspond to all the rows or columns.
5010 //The use of the list storage scheme, should give fast insertion
5011 //and fast sorts later.
5013  matrix_data_list(n_matrix);
5014  //Loop over the number of matrices and resize
5015  for(unsigned m=0;m<n_matrix;m++) {matrix_data_list[m].resize(ndof);}
5016 
5017  //Resize the residuals vectors
5018  for(unsigned v=0;v<n_vector;v++)
5019  {
5020  residuals[v] = new double[ndof];
5021  for (unsigned i = 0; i < ndof; i++)
5022  {
5023  residuals[v][i] = 0;
5024  }
5025  }
5026 
5027 #ifdef OOMPH_HAS_MPI
5028 
5029 
5030  // Storage for assembly time for elements
5031  double t_assemble_start=0.0;
5032 
5033  // Storage for assembly times
5034  if ((!doing_residuals)&&
5036  {
5037  Elemental_assembly_time.resize(n_elements);
5038  }
5039 
5040 #endif
5041 
5042  //------------Assemble and populate the lists-----------------------
5043  {
5044  //Allocate local storage for the element's contribution to the
5045  //residuals vectors and system matrices of the size of the maximum
5046  //number of dofs in any element.
5047  //This means that the stored is only allocated (and deleted) once
5048  Vector<Vector<double> > el_residuals(n_vector);
5049  Vector<DenseMatrix<double> > el_jacobian(n_matrix);
5050 
5051 
5052  //Pointer to a single list to be used during the assembly
5053  std::list<std::pair<unsigned,double> > *list_pt;
5054 
5055  //Loop over the all elements
5056  for(unsigned long e=el_lo;e<=el_hi;e++)
5057  {
5058 
5059 #ifdef OOMPH_HAS_MPI
5060  // Time it?
5061  if ((!doing_residuals)&&
5063  {
5064  t_assemble_start=TimingHelpers::timer();
5065  }
5066 #endif
5067 
5068  //Get the pointer to the element
5069  GeneralisedElement* elem_pt = mesh_pt()->element_pt(e);
5070 
5071 #ifdef OOMPH_HAS_MPI
5072  //Ignore halo elements
5073  if (!elem_pt->is_halo())
5074  {
5075 #endif
5076 
5077  //Find number of degrees of freedom in the element
5078  const unsigned nvar = assembly_handler_pt->ndof(elem_pt);
5079 
5080  //Resize the storage for the elemental jacobian and residuals
5081  for(unsigned v=0;v<n_vector;v++) {el_residuals[v].resize(nvar);}
5082  for(unsigned m=0;m<n_matrix;m++) {el_jacobian[m].resize(nvar);}
5083 
5084  //Now get the residuals and jacobian for the element
5085  assembly_handler_pt->
5086  get_all_vectors_and_matrices(elem_pt,el_residuals, el_jacobian);
5087 
5088  //---------------- Insert the values into the lists -----------
5089 
5090  //Loop over the first index of local variables
5091  for(unsigned i=0;i<nvar;i++)
5092  {
5093  //Get the local equation number
5094  unsigned eqn_number
5095  = assembly_handler_pt->eqn_number(elem_pt,i);
5096 
5097  //Add the contribution to the residuals
5098  for(unsigned v=0;v<n_vector;v++)
5099  {
5100  //Fill in the residuals vector
5101  residuals[v][eqn_number] += el_residuals[v][i];
5102  }
5103 
5104  //Now loop over the other index
5105  for(unsigned j=0;j<nvar;j++)
5106  {
5107  //Get the number of the unknown
5108  unsigned unknown = assembly_handler_pt->eqn_number(elem_pt,j);
5109 
5110  //Loop over the matrices
5111  for(unsigned m=0;m<n_matrix;m++)
5112  {
5113  //Get the value of the matrix at this point
5114  double value = el_jacobian[m](i,j);
5115  //Only add to theif it's non-zero
5116  if(std::fabs(value) > Numerical_zero_for_sparse_assembly)
5117  {
5118  //If it's compressed row storage, then our vector is indexed
5119  //by row (the equation number)
5120  if(compressed_row_flag)
5121  {
5122  //Find the list that corresponds to the desired row
5123  list_pt = &matrix_data_list[m][eqn_number];
5124  //Insert the data into the list, the first entry
5125  //in the pair is the unknown (column index),
5126  //the second is the value itself.
5127  list_pt->
5128  insert(list_pt->end(),std::make_pair(unknown,value));
5129  }
5130  //Otherwise it's compressed column storage, and our
5131  //vector is indexed by column (the unknown)
5132  else
5133  {
5134  //Find the list that corresponds to the desired column
5135  list_pt = &matrix_data_list[m][unknown];
5136  //Insert the data into the list, the first entry
5137  //in the pair is the equation number (row index),
5138  //the second is the value itself.
5139  list_pt->
5140  insert(list_pt->end(),std::make_pair(eqn_number,value));
5141  }
5142  }
5143  }
5144  }
5145  }
5146 
5147 #ifdef OOMPH_HAS_MPI
5148  } // endif halo element
5149 #endif
5150 
5151 
5152 #ifdef OOMPH_HAS_MPI
5153  // Time it?
5154  if ((!doing_residuals)&&
5156  {
5157  Elemental_assembly_time[e]=TimingHelpers::timer()-t_assemble_start;
5158  }
5159 #endif
5160 
5161  } //End of loop over the elements
5162 
5163  } //list_pt goes out of scope
5164 
5165 
5166 #ifdef OOMPH_HAS_MPI
5167 
5168  // Postprocess timing information and re-allocate distribution of
5169  // elements during subsequent assemblies.
5170  if ((!doing_residuals)&&
5173  {
5175  }
5176 
5177  // We have determined load balancing for current setup.
5178  // This can remain the same until assign_eqn_numbers() is called
5179  // again -- the flag is re-set to true there.
5180  if ((!doing_residuals)&&
5182  {
5184  }
5185 
5186 #endif
5187 
5188 
5189  //----Finally we need to convert the beautiful list storage scheme---
5190  //----------to the containers required by SuperLU--------------------
5191 
5192  //Loop over the number of matrices
5193  for(unsigned m=0;m<n_matrix;m++)
5194  {
5195  //Set the number of rows or columns
5196  row_or_column_start[m] = new int[ndof+1];
5197  //Counter for the total number of entries in the storage scheme
5198  unsigned long entry_count=0;
5199  //The first entry is 0
5200  row_or_column_start[m][0] = entry_count;
5201 
5202  // first we compute the number of non-zeros
5203  nnz[m] = 0;
5204  for(unsigned long i_global=0;i_global<ndof;i_global++)
5205  {
5206  nnz[m] += matrix_data_list[m][i_global].size();
5207  }
5208 
5209  // and then resize the storage
5210  column_or_row_index[m] = new int[nnz[m]];
5211  value[m] = new double[nnz[m]];
5212 
5213  //Now we merely loop over the number of rows or columns
5214  for(unsigned long i_global=0;i_global<ndof;i_global++)
5215  {
5216  //Start index for the present row is the number of entries so far
5217  row_or_column_start[m][i_global] = entry_count;
5218  //If there are no entries in the list then skip the loop
5219  if(matrix_data_list[m][i_global].empty()) {continue;}
5220 
5221  //Sort the list corresponding to this row or column by the
5222  //column or row index (first entry in the pair).
5223  //This might be inefficient, but we only have to do the sort ONCE
5224  //for each list. This is faster than using a map storage scheme, where
5225  //we are sorting for every insertion (although the map structure
5226  //is cleaner and more memory efficient)
5227  matrix_data_list[m][i_global].sort();
5228 
5229  //Set up an iterator for start of the list
5230  std::list<std::pair<unsigned,double> >::iterator it
5231  = matrix_data_list[m][i_global].begin();
5232 
5233  //Get the first row or column index in the list...
5234  unsigned current_index = it->first;
5235  //...and the corresponding value
5236  double current_value = it->second;
5237 
5238  //Loop over all the entries in the sorted list
5239  //Increase the iterator so that we start at the second entry
5240  for(++it;it!=matrix_data_list[m][i_global].end();++it)
5241  {
5242  //If the index has not changed, then we must add the contribution
5243  //of the present entry to the value.
5244  //Additionally check that the entry is non-zero
5245  if((it->first == current_index) &&
5246  (std::fabs(it->second) > Numerical_zero_for_sparse_assembly))
5247  {
5248  current_value += it->second;
5249  }
5250  //Otherwise, we have added all the contributions to the index
5251  //to current_value, so add it to the SuperLU data structure
5252  else
5253  {
5254  //Add the row or column index to the vector
5255  column_or_row_index[m][entry_count] = current_index;
5256  //Add the actual value to the vector
5257  value[m][entry_count] = current_value;
5258  //Increase the counter for the number of entries in each vector
5259  entry_count++;
5260 
5261  //Set the index and value to be those of the current entry in the
5262  //list
5263  current_index = it->first;
5264  current_value = it->second;
5265  }
5266  } //End of loop over all list entries for this global row or column
5267 
5268  //There are TWO special cases to consider.
5269  //If there is only one equation number in the list, then it
5270  //will NOT have been added. We test this case by comparing the
5271  //number of entries with those stored in row_or_column_start[i_global]
5272  //Otherwise
5273  //If the final entry in the list has the same index as the penultimate
5274  //entry, then it will NOT have been added to the SuperLU storage scheme
5275  //Check this by comparing the current_index with the final index
5276  //stored in the SuperLU scheme. If they are not the same, then
5277  //add the current_index and value.
5278 
5279  //If single equation number in list
5280  if((static_cast<int>(entry_count) == row_or_column_start[m][i_global])
5281  //If we have a single equation number, this will not be evaluated.
5282  //If we don't then we do the test to check that the final
5283  //entry is added
5284  ||(static_cast<int>(current_index) !=
5285  column_or_row_index[m][entry_count-1]))
5286  {
5287  //Add the row or column index to the vector
5288  column_or_row_index[m][entry_count] = current_index;
5289  //Add the actual value to the vector
5290  value[m][entry_count] = current_value;
5291  //Increase the counter for the number of entries in each vector
5292  entry_count++;
5293  }
5294 
5295  } //End of loop over the rows or columns of the entire matrix
5296 
5297  //Final entry in the row/column start vector
5298  row_or_column_start[m][ndof] = entry_count;
5299  } //End of loop over matrices
5300 
5302  {
5303  oomph_info << "Pausing at end of sparse assembly." << std::endl;
5304  pause("Check memory usage now.");
5305  }
5306 
5307 }
5308 
5309 
5310 
5311 //=====================================================================
5312 /// This is a (private) helper function that is used to assemble system
5313 /// matrices in compressed row or column format
5314 /// and compute residual vectors using vectors of pairs
5315 /// The default action is to assemble the jacobian matrix and
5316 /// residuals for the Newton method. The action can be
5317 /// overloaded at an elemental level by chaging the default
5318 /// behaviour of the function Element::get_all_vectors_and_matrices().
5319 /// column_or_row_index: Column [or row] index of given entry
5320 /// row_or_column_start: Index of first entry for given row [or column]
5321 /// value : Vector of nonzero entries
5322 /// residuals : Residual vector
5323 /// compressed_row_flag: Bool flag to indicate if storage format is
5324 /// compressed row [if false interpretation of
5325 /// arguments is as stated in square brackets].
5326 //=====================================================================
5328  Vector<int* > &column_or_row_index,
5329  Vector<int* > &row_or_column_start,
5330  Vector<double* > &value,
5331  Vector<unsigned> &nnz,
5332  Vector<double*> &residuals,
5333  bool compressed_row_flag)
5334 {
5335  //Total number of elements
5336  const unsigned long n_elements = mesh_pt()->nelement();
5337 
5338  // Default range of elements for distributed problems
5339  unsigned long el_lo=0;
5340  unsigned long el_hi=n_elements-1;
5341 
5342 #ifdef OOMPH_HAS_MPI
5343  // Otherwise just loop over a fraction of the elements
5344  // (This will either have been initialised in
5345  // Problem::set_default_first_and_last_element_for_assembly() or
5346  // will have been re-assigned during a previous assembly loop
5347  // Note that following the re-assignment only the entries
5348  // for the current processor are relevant.
5350  {
5351  el_lo=First_el_for_assembly[Communicator_pt->my_rank()];
5352  el_hi=Last_el_plus_one_for_assembly[Communicator_pt->my_rank()]-1;
5353  }
5354 #endif
5355 
5356  // number of local eqns
5357  unsigned ndof = this->ndof();
5358 
5359  //Find the number of vectors to be assembled
5360  const unsigned n_vector = residuals.size();
5361 
5362  //Find the number of matrices to be assembled
5363  const unsigned n_matrix = column_or_row_index.size();
5364 
5365  //Locally cache pointer to assembly handler
5367 
5368 #ifdef OOMPH_HAS_MPI
5369  bool doing_residuals=false;
5370  if (dynamic_cast<ParallelResidualsHandler*>(Assembly_handler_pt)!=0)
5371  {
5372  doing_residuals=true;
5373  }
5374 #endif
5375 
5376 //Error check dimensions
5377 #ifdef PARANOID
5378  if(row_or_column_start.size() != n_matrix)
5379  {
5380  std::ostringstream error_stream;
5381  error_stream
5382  << "Error: " << std::endl
5383  << "row_or_column_start.size() " << row_or_column_start.size()
5384  << " does not equal "
5385  << "column_or_row_index.size() "
5386  << column_or_row_index.size() << std::endl;
5387  throw OomphLibError(
5388  error_stream.str(),
5389  OOMPH_CURRENT_FUNCTION,
5390  OOMPH_EXCEPTION_LOCATION);
5391  }
5392 
5393  if(value.size() != n_matrix)
5394  {
5395  std::ostringstream error_stream;
5396  error_stream
5397  << "Error: "
5398  << std::endl
5399  << "value.size() " << value.size() << " does not equal "
5400  << "column_or_row_index.size() "
5401  << column_or_row_index.size() << std::endl<< std::endl
5402  << std::endl;
5403  throw OomphLibError(
5404  error_stream.str(),
5405  OOMPH_CURRENT_FUNCTION,
5406  OOMPH_EXCEPTION_LOCATION);
5407  }
5408 #endif
5409 
5410 
5411 // The idea behind this sparse assembly routine is to use a Vector of
5412 // Vectors of pairs for each complete matrix.
5413 // Each inner Vector stores pairs and holds the row (or column) index
5414 // and the value of the matrix entry.
5415 
5416 // Set up Vector of Vectors to store the entries of each matrix,
5417 // indexed by either the column or row.
5418  Vector<Vector< Vector<std::pair<unsigned,double> > > > matrix_data(n_matrix);
5419 
5420  //Loop over the number of matrices being assembled and resize
5421  //each Vector of Vectors to the number of rows or columns of the matrix
5422  for(unsigned m=0;m<n_matrix;m++) {matrix_data[m].resize(ndof);}
5423 
5424  //Resize the residuals vectors
5425  for(unsigned v=0;v<n_vector;v++)
5426  {
5427  residuals[v] = new double[ndof];
5428  for (unsigned i = 0; i < ndof; i++)
5429  {
5430  residuals[v][i] = 0;
5431  }
5432  }
5433 
5434 #ifdef OOMPH_HAS_MPI
5435 
5436  // Storage for assembly time for elements
5437  double t_assemble_start=0.0;
5438 
5439  // Storage for assembly times
5440  if ((!doing_residuals)&&
5442  {
5443  Elemental_assembly_time.resize(n_elements);
5444  }
5445 
5446 #endif
5447 
5448  //----------------Assemble and populate the vector storage scheme--------
5449  {
5450  //Allocate local storage for the element's contribution to the
5451  //residuals vectors and system matrices of the size of the maximum
5452  //number of dofs in any element
5453  //This means that the storage is only allocated (and deleted) once
5454  Vector<Vector<double> > el_residuals(n_vector);
5455  Vector<DenseMatrix<double> > el_jacobian(n_matrix);
5456 
5457  //Loop over the elements
5458  for(unsigned long e=el_lo;e<=el_hi;e++)
5459  {
5460 
5461 #ifdef OOMPH_HAS_MPI
5462  // Time it?
5463  if ((!doing_residuals)&&
5465  {
5466  t_assemble_start=TimingHelpers::timer();
5467  }
5468 #endif
5469 
5470  //Get the pointer to the element
5471  GeneralisedElement* elem_pt = mesh_pt()->element_pt(e);
5472 
5473 #ifdef OOMPH_HAS_MPI
5474  //Ignore halo elements
5475  if (!elem_pt->is_halo())
5476  {
5477 #endif
5478 
5479  //Find number of degrees of freedom in the element
5480  const unsigned nvar = assembly_handler_pt->ndof(elem_pt);
5481 
5482  //Resize the storage for elemental jacobian and residuals
5483  for(unsigned v=0;v<n_vector;v++) {el_residuals[v].resize(nvar);}
5484  for(unsigned m=0;m<n_matrix;m++) {el_jacobian[m].resize(nvar);}
5485 
5486  //Now get the residuals and jacobian for the element
5487  assembly_handler_pt->
5488  get_all_vectors_and_matrices(elem_pt,el_residuals, el_jacobian);
5489 
5490  //---------------Insert the values into the vectors--------------
5491 
5492  //Loop over the first index of local variables
5493  for(unsigned i=0;i<nvar;i++)
5494  {
5495 
5496  //Get the local equation number
5497  unsigned eqn_number
5498  = assembly_handler_pt->eqn_number(elem_pt,i);
5499 
5500  //Add the contribution to the residuals
5501  for(unsigned v=0;v<n_vector;v++)
5502  {
5503  //Fill in each residuals vector
5504  residuals[v][eqn_number] += el_residuals[v][i];
5505  }
5506 
5507  //Now loop over the other index
5508  for(unsigned j=0;j<nvar;j++)
5509  {
5510  //Get the number of the unknown
5511  unsigned unknown = assembly_handler_pt->eqn_number(elem_pt,j);
5512 
5513  //Loop over the matrices
5514  //If it's compressed row storage, then our vector of maps
5515  //is indexed by row (equation number)
5516  for(unsigned m=0;m<n_matrix;m++)
5517  {
5518  //Get the value of the matrix at this point
5519  double value = el_jacobian[m](i,j);
5520  //Only bother to add to the vector if it's non-zero
5521  if(std::fabs(value) > Numerical_zero_for_sparse_assembly)
5522  {
5523  //If it's compressed row storage, then our vector of maps
5524  //is indexed by row (equation number)
5525  if(compressed_row_flag)
5526  {
5527  //Find the correct position and add the data into the vectors
5528  const unsigned size = matrix_data[m][eqn_number].size();
5529  for(unsigned k=0; k<=size; k++)
5530  {
5531  if(k==size)
5532  {
5533  matrix_data[m][eqn_number].push_back(
5534  std::make_pair(unknown,value));
5535  break;
5536  }
5537  else if(matrix_data[m][eqn_number][k].first == unknown)
5538  {
5539  matrix_data[m][eqn_number][k].second += value;
5540  break;
5541  }
5542  }
5543  }
5544  //Otherwise it's compressed column storage and our vector is
5545  //indexed by column (the unknown)
5546  else
5547  {
5548  //Add the data into the vectors in the correct position
5549  const unsigned size = matrix_data[m][unknown].size();
5550  for(unsigned k=0; k<=size; k++)
5551  {
5552  if(k==size)
5553  {
5554  matrix_data[m][unknown].push_back(
5555  std::make_pair(eqn_number,value));
5556  break;
5557  }
5558  else if(matrix_data[m][unknown][k].first == eqn_number)
5559  {
5560  matrix_data[m][unknown][k].second += value;
5561  break;
5562  }
5563  }
5564  }
5565  }
5566  } //End of loop over matrices
5567  }
5568  }
5569 
5570 #ifdef OOMPH_HAS_MPI
5571  } // endif halo element
5572 #endif
5573 
5574 
5575 #ifdef OOMPH_HAS_MPI
5576  // Time it?
5577  if ((!doing_residuals)&&
5579  {
5580  Elemental_assembly_time[e]=TimingHelpers::timer()-t_assemble_start;
5581  }
5582 #endif
5583 
5584  } //End of loop over the elements
5585 
5586 
5587  } //End of vector assembly
5588 
5589 
5590 
5591 #ifdef OOMPH_HAS_MPI
5592 
5593  // Postprocess timing information and re-allocate distribution of
5594  // elements during subsequent assemblies.
5595  if ((!doing_residuals)&&
5598  {
5600  }
5601 
5602  // We have determined load balancing for current setup.
5603  // This can remain the same until assign_eqn_numbers() is called
5604  // again -- the flag is re-set to true there.
5605  if ((!doing_residuals)&&
5607  {
5609  }
5610 
5611 #endif
5612 
5613 
5614  //-----------Finally we need to convert this vector storage scheme
5615  //------------------------to the containers required by SuperLU
5616 
5617  //Loop over the number of matrices
5618  for(unsigned m=0;m<n_matrix;m++)
5619  {
5620  //Set the number of rows or columns
5621  row_or_column_start[m] = new int[ndof+1];
5622 
5623  // fill row_or_column_start and find the number of entries
5624  row_or_column_start[m][0] = 0;
5625  for(unsigned long i=0;i<ndof;i++)
5626  {
5627  row_or_column_start[m][i+1] = row_or_column_start[m][i]
5628  + matrix_data[m][i].size();
5629  }
5630  const unsigned entries = row_or_column_start[m][ndof];
5631 
5632  // resize vectors
5633  column_or_row_index[m] = new int[entries];
5634  value[m] = new double[entries];
5635  nnz[m] = entries;
5636 
5637  //Now we merely loop over the number of rows or columns
5638  for(unsigned long i_global=0;i_global<ndof;i_global++)
5639  {
5640  //If there are no entries in the vector then skip the rest of the loop
5641  if(matrix_data[m][i_global].empty()) {continue;}
5642 
5643  //Loop over all the entries in the vectors corresponding to the given
5644  //row or column. It will NOT be ordered
5645  unsigned p = 0;
5646  for(int j=row_or_column_start[m][i_global];
5647  j<row_or_column_start[m][i_global+1]; j++)
5648  {
5649  column_or_row_index[m][j] = matrix_data[m][i_global][p].first;
5650  value[m][j] = matrix_data[m][i_global][p].second;
5651  ++p;
5652  }
5653  }
5654  } //End of the loop over the matrices
5655 
5657  {
5658  oomph_info << "Pausing at end of sparse assembly." << std::endl;
5659  pause("Check memory usage now.");
5660  }
5661 }
5662 
5663 
5664 
5665 
5666 
5667 
5668 
5669 
5670 //=====================================================================
5671 /// This is a (private) helper function that is used to assemble system
5672 /// matrices in compressed row or column format
5673 /// and compute residual vectors using two vectors.
5674 /// The default action is to assemble the jacobian matrix and
5675 /// residuals for the Newton method. The action can be
5676 /// overloaded at an elemental level by chaging the default
5677 /// behaviour of the function Element::get_all_vectors_and_matrices().
5678 /// column_or_row_index: Column [or row] index of given entry
5679 /// row_or_column_start: Index of first entry for given row [or column]
5680 /// value : Vector of nonzero entries
5681 /// residuals : Residual vector
5682 /// compressed_row_flag: Bool flag to indicate if storage format is
5683 /// compressed row [if false interpretation of
5684 /// arguments is as stated in square brackets].
5685 //=====================================================================
5687  Vector<int* > &column_or_row_index,
5688  Vector<int* > &row_or_column_start,
5689  Vector<double* > &value,
5690  Vector<unsigned> &nnz,
5691  Vector<double* > &residuals,
5692  bool compressed_row_flag)
5693 {
5694  //Total number of elements
5695  const unsigned long n_elements = mesh_pt()->nelement();
5696 
5697  // Default range of elements for distributed problems
5698  unsigned long el_lo=0;
5699  unsigned long el_hi=n_elements-1;
5700 
5701 
5702 #ifdef OOMPH_HAS_MPI
5703  // Otherwise just loop over a fraction of the elements
5704  // (This will either have been initialised in
5705  // Problem::set_default_first_and_last_element_for_assembly() or
5706  // will have been re-assigned during a previous assembly loop
5707  // Note that following the re-assignment only the entries
5708  // for the current processor are relevant.
5710  {
5711  el_lo=First_el_for_assembly[Communicator_pt->my_rank()];
5712  el_hi=Last_el_plus_one_for_assembly[Communicator_pt->my_rank()]-1;
5713  }
5714 #endif
5715 
5716  // number of local eqns
5717  unsigned ndof = this->ndof();
5718 
5719  //Find the number of vectors to be assembled
5720  const unsigned n_vector = residuals.size();
5721 
5722  //Find the number of matrices to be assembled
5723  const unsigned n_matrix = column_or_row_index.size();
5724 
5725  //Locally cache pointer to assembly handler
5727 
5728 #ifdef OOMPH_HAS_MPI
5729  bool doing_residuals=false;
5730  if (dynamic_cast<ParallelResidualsHandler*>(Assembly_handler_pt)!=0)
5731  {
5732  doing_residuals=true;
5733  }
5734 #endif
5735 
5736 //Error check dimensions
5737 #ifdef PARANOID
5738  if(row_or_column_start.size() != n_matrix)
5739  {
5740  std::ostringstream error_stream;
5741  error_stream
5742  << "Error: " << std::endl
5743  << "row_or_column_start.size() " << row_or_column_start.size()
5744  << " does not equal "
5745  << "column_or_row_index.size() "
5746  << column_or_row_index.size() << std::endl;
5747  throw OomphLibError(
5748  error_stream.str(),
5749  OOMPH_CURRENT_FUNCTION,
5750  OOMPH_EXCEPTION_LOCATION);
5751  }
5752 
5753  if(value.size() != n_matrix)
5754  {
5755  std::ostringstream error_stream;
5756  error_stream
5757  << "Error: "
5758  << std::endl
5759  << "value.size() " << value.size() << " does not equal "
5760  << "column_or_row_index.size() "
5761  << column_or_row_index.size() << std::endl<< std::endl
5762  << std::endl;
5763  throw OomphLibError(
5764  error_stream.str(),
5765  OOMPH_CURRENT_FUNCTION,
5766  OOMPH_EXCEPTION_LOCATION);
5767  }
5768 #endif
5769 
5770 // The idea behind this sparse assembly routine is to use Vectors of
5771 // Vectors for the entries in each complete matrix. And a second
5772 // Vector of Vectors stores the global row (or column) indeces. This
5773 // will not have the memory overheads associated with the methods using
5774 // lists or maps, but insertion will be more costly.
5775 
5776 // Set up two vector of vectors to store the entries of each matrix,
5777 // indexed by either the column or row. The entries of the vector for
5778 // each matrix correspond to all the rows or columns of that matrix.
5779  Vector<Vector<Vector<unsigned> > > matrix_row_or_col_indices(n_matrix);
5780  Vector<Vector<Vector<double> > > matrix_values(n_matrix);
5781 
5782  //Loop over the number of matrices being assembled and resize
5783  //each vector of vectors to the number of rows or columns of the matrix
5784  for(unsigned m=0;m<n_matrix;m++)
5785  {
5786  matrix_row_or_col_indices[m].resize(ndof);
5787  matrix_values[m].resize(ndof);
5788  }
5789 
5790  //Resize the residuals vectors
5791  for(unsigned v=0;v<n_vector;v++)
5792  {
5793  residuals[v] = new double[ndof];
5794  for (unsigned i = 0; i < ndof; i++)
5795  {
5796  residuals[v][i] = 0;
5797  }
5798  }
5799 
5800 #ifdef OOMPH_HAS_MPI
5801 
5802  // Storage for assembly time for elements
5803  double t_assemble_start=0.0;
5804 
5805  // Storage for assembly times
5806  if ((!doing_residuals)&&
5808  {
5809  Elemental_assembly_time.resize(n_elements);
5810  }
5811 
5812 #endif
5813 
5814 
5815  //----------------Assemble and populate the vector storage scheme-------
5816  {
5817  //Allocate local storage for the element's contribution to the
5818  //residuals vectors and system matrices of the size of the maximum
5819  //number of dofs in any element
5820  //This means that the storage will only be allocated (and deleted) once
5821  Vector<Vector<double> > el_residuals(n_vector);
5822  Vector<DenseMatrix<double> > el_jacobian(n_matrix);
5823 
5824  //Loop over the elements
5825  for(unsigned long e=el_lo;e<=el_hi;e++)
5826  {
5827 
5828 #ifdef OOMPH_HAS_MPI
5829  // Time it?
5830  if ((!doing_residuals)&&
5832  {
5833  t_assemble_start=TimingHelpers::timer();
5834  }
5835 #endif
5836 
5837  //Get the pointer to the element
5838  GeneralisedElement* elem_pt = mesh_pt()->element_pt(e);
5839 
5840 #ifdef OOMPH_HAS_MPI
5841  //Ignore halo elements
5842  if (!elem_pt->is_halo())
5843  {
5844 #endif
5845 
5846  //Find number of degrees of freedom in the element
5847  const unsigned nvar = assembly_handler_pt->ndof(elem_pt);
5848 
5849  //Resize the storage for elemental jacobian and residuals
5850  for(unsigned v=0;v<n_vector;v++) {el_residuals[v].resize(nvar);}
5851  for(unsigned m=0;m<n_matrix;m++) {el_jacobian[m].resize(nvar);}
5852 
5853  //Now get the residuals and jacobian for the element
5854  assembly_handler_pt->
5855  get_all_vectors_and_matrices(elem_pt,el_residuals, el_jacobian);
5856 
5857  //---------------Insert the values into the vectors--------------
5858 
5859  //Loop over the first index of local variables
5860  for(unsigned i=0;i<nvar;i++)
5861  {
5862  //Get the local equation number
5863  unsigned eqn_number
5864  = assembly_handler_pt->eqn_number(elem_pt,i);
5865 
5866  //Add the contribution to the residuals
5867  for(unsigned v=0;v<n_vector;v++)
5868  {
5869  //Fill in each residuals vector
5870  residuals[v][eqn_number] += el_residuals[v][i];
5871  }
5872 
5873  //Now loop over the other index
5874  for(unsigned j=0;j<nvar;j++)
5875  {
5876  //Get the number of the unknown
5877  unsigned unknown = assembly_handler_pt->eqn_number(elem_pt,j);
5878 
5879  //Loop over the matrices
5880  //If it's compressed row storage, then our vector of maps
5881  //is indexed by row (equation number)
5882  for(unsigned m=0;m<n_matrix;m++)
5883  {
5884  //Get the value of the matrix at this point
5885  double value = el_jacobian[m](i,j);
5886  //Only bother to add to the vector if it's non-zero
5887  if(std::fabs(value) > Numerical_zero_for_sparse_assembly)
5888  {
5889  //If it's compressed row storage, then our vector of maps
5890  //is indexed by row (equation number)
5891  if(compressed_row_flag)
5892  {
5893  //Find the correct position and add the data into the vectors
5894  const unsigned size =
5895  matrix_row_or_col_indices[m][eqn_number].size();
5896 
5897  for(unsigned k=0; k<=size; k++)
5898  {
5899  if(k==size)
5900  {
5901  matrix_row_or_col_indices[m][eqn_number].
5902  push_back(unknown);
5903  matrix_values[m][eqn_number].push_back(value);
5904  break;
5905  }
5906  else if(matrix_row_or_col_indices[m][eqn_number][k] ==
5907  unknown)
5908  {
5909  matrix_values[m][eqn_number][k] += value;
5910  break;
5911  }
5912  }
5913  }
5914  //Otherwise it's compressed column storage and our vector is
5915  //indexed by column (the unknown)
5916  else
5917  {
5918  //Add the data into the vectors in the correct position
5919  const unsigned size =
5920  matrix_row_or_col_indices[m][unknown].size();
5921  for (unsigned k=0; k<=size; k++)
5922  {
5923  if (k==size)
5924  {
5925  matrix_row_or_col_indices[m][unknown].
5926  push_back(eqn_number);
5927  matrix_values[m][unknown].push_back(value);
5928  break;
5929  }
5930  else if (matrix_row_or_col_indices[m][unknown][k] ==
5931  eqn_number)
5932  {
5933  matrix_values[m][unknown][k] += value;
5934  break;
5935  }
5936  }
5937  }
5938  }
5939  } //End of loop over matrices
5940  }
5941  }
5942 
5943 #ifdef OOMPH_HAS_MPI
5944  } // endif halo element
5945 #endif
5946 
5947 
5948 #ifdef OOMPH_HAS_MPI
5949  // Time it?
5950  if ((!doing_residuals)&&
5952  {
5953  Elemental_assembly_time[e]=TimingHelpers::timer()-t_assemble_start;
5954  }
5955 #endif
5956 
5957  } //End of loop over the elements
5958 
5959  } //End of vector assembly
5960 
5961 
5962 
5963 #ifdef OOMPH_HAS_MPI
5964 
5965  // Postprocess timing information and re-allocate distribution of
5966  // elements during subsequent assemblies.
5967  if ((!doing_residuals)&&
5970  {
5972  }
5973 
5974  // We have determined load balancing for current setup.
5975  // This can remain the same until assign_eqn_numbers() is called
5976  // again -- the flag is re-set to true there.
5977  if ((!doing_residuals)&&
5979  {
5981  }
5982 
5983 #endif
5984 
5985  //-----------Finally we need to convert this lousy vector storage scheme
5986  //------------------------to the containers required by SuperLU
5987 
5988  //Loop over the number of matrices
5989  for(unsigned m=0;m<n_matrix;m++)
5990  {
5991  //Set the number of rows or columns
5992  row_or_column_start[m] = new int[ndof+1];
5993 
5994  // fill row_or_column_start and find the number of entries
5995  row_or_column_start[m][0] = 0;
5996  for (unsigned long i=0;i<ndof;i++)
5997  {
5998  row_or_column_start[m][i+1] = row_or_column_start[m][i]
5999  + matrix_values[m][i].size();
6000  }
6001  const unsigned entries = row_or_column_start[m][ndof];
6002 
6003  // resize vectors
6004  column_or_row_index[m] = new int[entries];
6005  value[m] = new double[entries];
6006  nnz[m] = entries;
6007 
6008  //Now we merely loop over the number of rows or columns
6009  for(unsigned long i_global=0;i_global<ndof;i_global++)
6010  {
6011  //If there are no entries in the vector then skip the rest of the loop
6012  if(matrix_values[m][i_global].empty()) {continue;}
6013 
6014  //Loop over all the entries in the vectors corresponding to the given
6015  //row or column. It will NOT be ordered
6016  unsigned p = 0;
6017  for(int j=row_or_column_start[m][i_global];
6018  j<row_or_column_start[m][i_global+1]; j++)
6019  {
6020  column_or_row_index[m][j] = matrix_row_or_col_indices[m][i_global][p];
6021  value[m][j] = matrix_values[m][i_global][p];
6022  ++p;
6023  }
6024  }
6025  } //End of the loop over the matrices
6026 
6028  {
6029  oomph_info << "Pausing at end of sparse assembly." << std::endl;
6030  pause("Check memory usage now.");
6031  }
6032 
6033 }
6034 
6035 
6036 //=====================================================================
6037 /// This is a (private) helper function that is used to assemble system
6038 /// matrices in compressed row or column format
6039 /// and compute residual vectors using two vectors.
6040 /// The default action is to assemble the jacobian matrix and
6041 /// residuals for the Newton method. The action can be
6042 /// overloaded at an elemental level by chaging the default
6043 /// behaviour of the function Element::get_all_vectors_and_matrices().
6044 /// column_or_row_index: Column [or row] index of given entry
6045 /// row_or_column_start: Index of first entry for given row [or column]
6046 /// value : Vector of nonzero entries
6047 /// residuals : Residual vector
6048 /// compressed_row_flag: Bool flag to indicate if storage format is
6049 /// compressed row [if false interpretation of
6050 /// arguments is as stated in square brackets].
6051 //=====================================================================
6053  Vector<int* > &column_or_row_index,
6054  Vector<int* > &row_or_column_start,
6055  Vector<double* > &value,
6056  Vector<unsigned> &nnz,
6057  Vector<double* > &residuals,
6058  bool compressed_row_flag)
6059 {
6060 
6061  //Total number of elements
6062  const unsigned long n_elements = mesh_pt()->nelement();
6063 
6064  // Default range of elements for distributed problems
6065  unsigned long el_lo=0;
6066  unsigned long el_hi=n_elements-1;
6067 
6068 
6069 #ifdef OOMPH_HAS_MPI
6070  // Otherwise just loop over a fraction of the elements
6071  // (This will either have been initialised in
6072  // Problem::set_default_first_and_last_element_for_assembly() or
6073  // will have been re-assigned during a previous assembly loop
6074  // Note that following the re-assignment only the entries
6075  // for the current processor are relevant.
6077  {
6078  el_lo=First_el_for_assembly[Communicator_pt->my_rank()];
6079  el_hi=Last_el_plus_one_for_assembly[Communicator_pt->my_rank()]-1;
6080  }
6081 #endif
6082 
6083  // number of local eqns
6084  unsigned ndof = this->ndof();
6085 
6086  //Find the number of vectors to be assembled
6087  const unsigned n_vector = residuals.size();
6088 
6089  //Find the number of matrices to be assembled
6090  const unsigned n_matrix = column_or_row_index.size();
6091 
6092  //Locally cache pointer to assembly handler
6094 
6095 #ifdef OOMPH_HAS_MPI
6096  bool doing_residuals=false;
6097  if (dynamic_cast<ParallelResidualsHandler*>(Assembly_handler_pt)!=0)
6098  {
6099  doing_residuals=true;
6100  }
6101 #endif
6102 
6103 //Error check dimensions
6104 #ifdef PARANOID
6105  if(row_or_column_start.size() != n_matrix)
6106  {
6107  std::ostringstream error_stream;
6108  error_stream
6109  << "Error: " << std::endl
6110  << "row_or_column_start.size() " << row_or_column_start.size()
6111  << " does not equal "
6112  << "column_or_row_index.size() "
6113  << column_or_row_index.size() << std::endl;
6114  throw OomphLibError(
6115  error_stream.str(),
6116  OOMPH_CURRENT_FUNCTION,
6117  OOMPH_EXCEPTION_LOCATION);
6118  }
6119 
6120  if(value.size() != n_matrix)
6121  {
6122  std::ostringstream error_stream;
6123  error_stream
6124  << "Error: "
6125  << std::endl
6126  << "value.size() " << value.size() << " does not equal "
6127  << "column_or_row_index.size() "
6128  << column_or_row_index.size() << std::endl<< std::endl
6129  << std::endl;
6130  throw OomphLibError(
6131  error_stream.str(),
6132  OOMPH_CURRENT_FUNCTION,
6133  OOMPH_EXCEPTION_LOCATION);
6134  }
6135 #endif
6136 
6137 // The idea behind this sparse assembly routine is to use Vectors of
6138 // Vectors for the entries in each complete matrix. And a second
6139 // Vector of Vectors stores the global row (or column) indeces. This
6140 // will not have the memory overheads associated with the methods using
6141 // lists or maps, but insertion will be more costly.
6142 
6143 // Set up two vector of vectors to store the entries of each matrix,
6144 // indexed by either the column or row. The entries of the vector for
6145 // each matrix correspond to all the rows or columns of that matrix.
6146  Vector<unsigned**> matrix_row_or_col_indices(n_matrix);
6147  Vector<double**> matrix_values(n_matrix);
6148 
6149  //Loop over the number of matrices being assembled and resize
6150  //each vector of vectors to the number of rows or columns of the matrix
6151  for(unsigned m=0;m<n_matrix;m++)
6152  {
6153  matrix_row_or_col_indices[m] = new unsigned*[ndof];
6154  matrix_values[m] = new double*[ndof];
6155  }
6156 
6157  //Resize the residuals vectors
6158  for(unsigned v=0;v<n_vector;v++)
6159  {
6160  residuals[v] = new double[ndof];
6161  for (unsigned i = 0; i < ndof; i++)
6162  {
6163  residuals[v][i] = 0;
6164  }
6165  }
6166 
6167 #ifdef OOMPH_HAS_MPI
6168 
6169  // Storage for assembly time for elements
6170  double t_assemble_start=0.0;
6171 
6172  // Storage for assembly times
6173  if ((!doing_residuals)&&
6175  {
6176  Elemental_assembly_time.resize(n_elements);
6177  }
6178 
6179 #endif
6180 
6181  // number of coefficients in each row
6182  Vector<Vector<unsigned> > ncoef(n_matrix);
6183  for (unsigned m = 0; m < n_matrix; m++)
6184  {
6185  ncoef[m].resize(ndof,0);
6186  }
6187 
6189  {
6191  for (unsigned m = 0; m < n_matrix; m++)
6192  {
6194  }
6195  }
6196 
6197  //----------------Assemble and populate the vector storage scheme-------
6198  {
6199  //Allocate local storage for the element's contribution to the
6200  //residuals vectors and system matrices of the size of the maximum
6201  //number of dofs in any element
6202  //This means that the storage will only be allocated (and deleted) once
6203  Vector<Vector<double> > el_residuals(n_vector);
6204  Vector<DenseMatrix<double> > el_jacobian(n_matrix);
6205 
6206  //Loop over the elements
6207  for(unsigned long e=el_lo;e<=el_hi;e++)
6208  {
6209 
6210 #ifdef OOMPH_HAS_MPI
6211  // Time it?
6212  if ((!doing_residuals)&&
6214  {
6215  t_assemble_start=TimingHelpers::timer();
6216  }
6217 #endif
6218 
6219  //Get the pointer to the element
6220  GeneralisedElement* elem_pt = mesh_pt()->element_pt(e);
6221 
6222 #ifdef OOMPH_HAS_MPI
6223  //Ignore halo elements
6224  if (!elem_pt->is_halo())
6225  {
6226 #endif
6227 
6228  //Find number of degrees of freedom in the element
6229  const unsigned nvar = assembly_handler_pt->ndof(elem_pt);
6230 
6231  //Resize the storage for elemental jacobian and residuals
6232  for(unsigned v=0;v<n_vector;v++) {el_residuals[v].resize(nvar);}
6233  for(unsigned m=0;m<n_matrix;m++) {el_jacobian[m].resize(nvar);}
6234 
6235  //Now get the residuals and jacobian for the element
6236  assembly_handler_pt->
6237  get_all_vectors_and_matrices(elem_pt,el_residuals, el_jacobian);
6238 
6239  //---------------Insert the values into the vectors--------------
6240 
6241  //Loop over the first index of local variables
6242  for(unsigned i=0;i<nvar;i++)
6243  {
6244  //Get the local equation number
6245  unsigned eqn_number
6246  = assembly_handler_pt->eqn_number(elem_pt,i);
6247 
6248  //Add the contribution to the residuals
6249  for(unsigned v=0;v<n_vector;v++)
6250  {
6251  //Fill in each residuals vector
6252  residuals[v][eqn_number] += el_residuals[v][i];
6253  }
6254 
6255  //Now loop over the other index
6256  for(unsigned j=0;j<nvar;j++)
6257  {
6258  //Get the number of the unknown
6259  unsigned unknown = assembly_handler_pt->eqn_number(elem_pt,j);
6260 
6261  //Loop over the matrices
6262  //If it's compressed row storage, then our vector of maps
6263  //is indexed by row (equation number)
6264  for(unsigned m=0;m<n_matrix;m++)
6265  {
6266  //Get the value of the matrix at this point
6267  double value = el_jacobian[m](i,j);
6268  //Only bother to add to the vector if it's non-zero
6269  if(std::fabs(value) > Numerical_zero_for_sparse_assembly)
6270  {
6271  // number of entrys in this row
6272  const unsigned size = ncoef[m][eqn_number];
6273 
6274  // if no data has been allocated for this row then allocate
6275  if (size == 0)
6276  {
6277  // do we have previous allocation data
6279  [m][eqn_number] != 0)
6280  {
6281  matrix_row_or_col_indices[m][eqn_number] =
6282  new unsigned
6284  [eqn_number]];
6285  matrix_values[m][eqn_number] =
6286  new double
6287  [Sparse_assemble_with_arrays_previous_allocation[m]
6288  [eqn_number]];
6289  }
6290  else
6291  {
6292  matrix_row_or_col_indices[m][eqn_number] =
6293  new unsigned
6295  matrix_values[m][eqn_number] =
6296  new double
6299  [eqn_number]=
6301  }
6302  }
6303 
6304  //If it's compressed row storage, then our vector of maps
6305  //is indexed by row (equation number)
6306  if(compressed_row_flag)
6307  {
6308  // next add the data
6309  for(unsigned k=0; k<=size; k++)
6310  {
6311  if(k==size)
6312  {
6313  // do we need to allocate more storage
6315  [m][eqn_number] == ncoef[m][eqn_number])
6316  {
6317  unsigned new_allocation = ncoef[m][eqn_number]+
6319  double* new_values = new double[new_allocation];
6320  unsigned* new_indices = new unsigned[new_allocation];
6321  for (unsigned c = 0; c < ncoef[m][eqn_number]; c++)
6322  {
6323  new_values[c] = matrix_values[m][eqn_number][c];
6324  new_indices[c] =
6325  matrix_row_or_col_indices[m][eqn_number][c];
6326  }
6327  delete[] matrix_values[m][eqn_number];
6328  delete[] matrix_row_or_col_indices[m][eqn_number];
6329  matrix_values[m][eqn_number]=new_values;
6330  matrix_row_or_col_indices[m][eqn_number]=new_indices;
6332  [m][eqn_number] = new_allocation;
6333  }
6334  // and now add the data
6335  unsigned entry = ncoef[m][eqn_number];
6336  ncoef[m][eqn_number]++;
6337  matrix_row_or_col_indices[m][eqn_number][entry] = unknown;
6338  matrix_values[m][eqn_number][entry] = value;
6339  break;
6340  }
6341  else if(matrix_row_or_col_indices[m][eqn_number][k] ==
6342  unknown)
6343  {
6344  matrix_values[m][eqn_number][k] += value;
6345  break;
6346  }
6347  }
6348  }
6349  //Otherwise it's compressed column storage and our vector is
6350  //indexed by column (the unknown)
6351  else
6352  {
6353  //Add the data into the vectors in the correct position
6354  for(unsigned k=0; k<=size; k++)
6355  {
6356  if(k==size)
6357  {
6358  // do we need to allocate more storage
6360  [m][unknown] == ncoef[m][unknown])
6361  {
6362  unsigned new_allocation = ncoef[m][unknown]+
6364  double* new_values = new double[new_allocation];
6365  unsigned* new_indices = new unsigned[new_allocation];
6366  for (unsigned c = 0; c < ncoef[m][unknown]; c++)
6367  {
6368  new_values[c] = matrix_values[m][unknown][c];
6369  new_indices[c] =
6370  matrix_row_or_col_indices[m][unknown][c];
6371  }
6372  delete[] matrix_values[m][unknown];
6373  delete[] matrix_row_or_col_indices[m][unknown];
6375  [m][unknown] = new_allocation;
6376  }
6377  // and now add the data
6378  unsigned entry = ncoef[m][unknown];
6379  ncoef[m][unknown]++;
6380  matrix_row_or_col_indices[m][unknown][entry] = eqn_number;
6381  matrix_values[m][unknown][entry] = value;
6382  break;
6383  }
6384  else if(matrix_row_or_col_indices[m][unknown][k] ==
6385  eqn_number)
6386  {
6387  matrix_values[m][unknown][k] += value;
6388  break;
6389  }
6390  }
6391  }
6392  }
6393  } //End of loop over matrices
6394  }
6395  }
6396 
6397 #ifdef OOMPH_HAS_MPI
6398  } // endif halo element
6399 #endif
6400 
6401 
6402 #ifdef OOMPH_HAS_MPI
6403  // Time it?
6404  if ((!doing_residuals)&&
6406  {
6407  Elemental_assembly_time[e]=TimingHelpers::timer()-t_assemble_start;
6408  }
6409 #endif
6410 
6411  } //End of loop over the elements
6412 
6413  } //End of vector assembly
6414 
6415 
6416 
6417 #ifdef OOMPH_HAS_MPI
6418 
6419  // Postprocess timing information and re-allocate distribution of
6420  // elements during subsequent assemblies.
6421  if ((!doing_residuals)&&
6424  {
6426  }
6427 
6428  // We have determined load balancing for current setup.
6429  // This can remain the same until assign_eqn_numbers() is called
6430  // again -- the flag is re-set to true there.
6431  if ((!doing_residuals)&&
6433  {
6435  }
6436 
6437 #endif
6438 
6439  //-----------Finally we need to convert this lousy vector storage scheme
6440  //------------------------to the containers required by SuperLU
6441 
6442  //Loop over the number of matrices
6443  for(unsigned m=0;m<n_matrix;m++)
6444  {
6445  //Set the number of rows or columns
6446  row_or_column_start[m] = new int[ndof+1];
6447 
6448  // fill row_or_column_start and find the number of entries
6449  row_or_column_start[m][0] = 0;
6450  for (unsigned long i=0;i<ndof;i++)
6451  {
6452  row_or_column_start[m][i+1] = row_or_column_start[m][i]
6453  + ncoef[m][i];
6455  }
6456  const unsigned entries = row_or_column_start[m][ndof];
6457 
6458  // resize vectors
6459  column_or_row_index[m] = new int[entries];
6460  value[m] = new double[entries];
6461  nnz[m] = entries;
6462 
6463  //Now we merely loop over the number of rows or columns
6464  for(unsigned long i_global=0;i_global<ndof;i_global++)
6465  {
6466  //If there are no entries in the vector then skip the rest of the loop
6467  if(ncoef[m][i_global]==0) {continue;}
6468 
6469  //Loop over all the entries in the vectors corresponding to the given
6470