mesh.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: 1282 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-01-16 08:27:53 +0000 (Mon, 16 Jan 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 //Non-inline member functions for general mesh classes
31 
32 #ifdef OOMPH_HAS_MPI
33 #include "mpi.h"
34 #endif
35 
36 #include<algorithm>
37 #include<limits.h>
38 #include <typeinfo>
39 
40 
41 //oomph-lib headers
42 #include "oomph_utilities.h"
43 #include "mesh.h"
44 #include "problem.h"
45 #include "elastic_problems.h"
46 #include "refineable_mesh.h"
47 #include "triangle_mesh.h"
48 #include "shape.h"
49 
50 namespace oomph
51 {
52 
53 //======================================================
54 /// The Steady Timestepper
55 //======================================================
57 
58 
59 //=======================================================================
60 /// Static boolean flag to control warning about mesh level timesteppers
61 //=======================================================================
63 
64 //=======================================================================
65 /// Merge meshes.
66 /// Note: This simply merges the meshes' elements and nodes (ignoring
67 /// duplicates; no boundary information etc. is created).
68 //=======================================================================
69  void Mesh::merge_meshes(const Vector<Mesh*>& sub_mesh_pt)
70  {
71  // No boundary lookup scheme is set up for the combined mesh
73 
74  //Number of submeshes
75  unsigned nsub_mesh=sub_mesh_pt.size();
76 
77  // Initialise element, node and boundary counters for global mesh
78  unsigned long n_element=0;
79  unsigned long n_node=0;
80  unsigned n_bound=0;
81 
82  // Loop over submeshes and get total number of elements, nodes and
83  // boundaries
84  for(unsigned imesh=0;imesh<nsub_mesh;imesh++)
85  {
86  n_element += sub_mesh_pt[imesh]->nelement();
87  n_node += sub_mesh_pt[imesh]->nnode();
88  n_bound += sub_mesh_pt[imesh]->nboundary();
89  }
90 
91  // Reserve storage for element and node pointers
92  Element_pt.clear();
93  Element_pt.reserve(n_element);
94  Node_pt.clear();
95  Node_pt.reserve(n_node);
96 
97  //Resize vector of vectors of nodes
98  Boundary_node_pt.clear();
99  Boundary_node_pt.resize(n_bound);
100 
101  // Sets of pointers to elements and nodes (to exlude duplicates -- they
102  // shouldn't occur anyway but if they do, they must only be added
103  // once in the global mesh to avoid trouble in the timestepping)
104  std::set<GeneralisedElement*> element_set_pt;
105  std::set<Node*> node_set_pt;
106 
107  //Counter for total number of boundaries in all the submeshes
108  unsigned ibound_global=0;
109  //Loop over the number of submeshes
110  for(unsigned imesh=0;imesh<nsub_mesh;imesh++)
111  {
112  //Loop over the elements of the submesh and add to vector
113  //duplicates are ignored
114  unsigned nel_before=0;
115  unsigned long n_element=sub_mesh_pt[imesh]->nelement();
116  for (unsigned long e=0;e<n_element;e++)
117  {
118  GeneralisedElement* el_pt=sub_mesh_pt[imesh]->element_pt(e);
119  element_set_pt.insert(el_pt);
120  // Was it a duplicate?
121  unsigned nel_now=element_set_pt.size();
122  if (nel_now==nel_before)
123  {
124  std::ostringstream warning_stream;
125  warning_stream <<"WARNING: " << std::endl
126  <<"Element " << e << " in submesh " << imesh
127  <<" is a duplicate \n and was ignored when assembling"
128  <<" combined mesh." << std::endl;
129  OomphLibWarning(warning_stream.str(),
130  "Mesh::Mesh(const Vector<Mesh*>&)",
131  OOMPH_EXCEPTION_LOCATION);
132  }
133  else
134  {
135  Element_pt.push_back(el_pt);
136  }
137  nel_before=nel_now;
138  }
139 
140  //Loop over the nodes of the submesh and add to vector
141  //duplicates are ignored
142  unsigned nnod_before=0;
143  unsigned long n_node=sub_mesh_pt[imesh]->nnode();
144  for (unsigned long n=0;n<n_node;n++)
145  {
146  Node* nod_pt=sub_mesh_pt[imesh]->node_pt(n);
147  node_set_pt.insert(nod_pt);
148  // Was it a duplicate?
149  unsigned nnod_now=node_set_pt.size();
150  if (nnod_now==nnod_before)
151  {
152  std::ostringstream warning_stream;
153  warning_stream<<"WARNING: " << std::endl
154  <<"Node " << n << " in submesh " << imesh
155  <<" is a duplicate \n and was ignored when assembling "
156  << "combined mesh." << std::endl;
157  OomphLibWarning(warning_stream.str(),
158  "Mesh::Mesh(const Vector<Mesh*>&)",
159  OOMPH_EXCEPTION_LOCATION);
160  }
161  else
162  {
163  Node_pt.push_back(nod_pt);
164  }
165  nnod_before=nnod_now;
166  }
167 
168  //Loop over the boundaries of the submesh
169  unsigned n_bound=sub_mesh_pt[imesh]->nboundary();
170  for (unsigned ibound=0;ibound<n_bound;ibound++)
171  {
172  //Loop over the number of nodes on the boundary and add to the
173  //global vector
174  unsigned long n_bound_node=sub_mesh_pt[imesh]->nboundary_node(ibound);
175  for (unsigned long n=0;n<n_bound_node;n++)
176  {
177  Boundary_node_pt[ibound_global].push_back(
178  sub_mesh_pt[imesh]->boundary_node_pt(ibound,n));
179  }
180  //Increase the number of the global boundary counter
181  ibound_global++;
182  }
183  } //End of loop over submeshes
184 
185  }
186 
187 
188 //========================================================
189 /// Remove the information about nodes stored on the
190 /// b-th boundary of the mesh
191 //========================================================
192 void Mesh::remove_boundary_nodes(const unsigned &b)
193 {
194  //Loop over all the nodes on the boundary and call
195  //their remove_from_boundary function
196  unsigned n_boundary_node = Boundary_node_pt[b].size();
197  for(unsigned n=0;n<n_boundary_node;n++)
198  {
200  }
201  //Clear the storage
202  Boundary_node_pt[b].clear();
203 }
204 
205 //=================================================================
206 /// Remove all information about mesh boundaries
207 //================================================================
209 {
210  //Loop over each boundary call remove_boundary_nodes
211  unsigned n_bound = Boundary_node_pt.size();
212  for(unsigned b=0;b<n_bound;b++) {remove_boundary_nodes(b);}
213  //Clear the storage
214  Boundary_node_pt.clear();
215 }
216 
217 //============================================================
218 /// Remove the node node_pt from the b-th boundary of the mesh
219 /// This function also removes the information from the Node
220 /// itself
221 //===========================================================
222 void Mesh::remove_boundary_node(const unsigned &b, Node* const &node_pt)
223 {
224  //Find the location of the node in the boundary
226  std::find(Boundary_node_pt[b].begin(),
227  Boundary_node_pt[b].end(),
228  node_pt);
229  //If the node is on this boundary
230  if(it!=Boundary_node_pt[b].end())
231  {
232  //Remove the node from the mesh's list of boundary nodes
233  Boundary_node_pt[b].erase(it);
234  //Now remove the node's boundary information
235  node_pt->remove_from_boundary(b);
236  }
237  //If not do nothing
238 }
239 
240 
241 //========================================================
242 /// Add the node node_pt to the b-th boundary of the mesh
243 /// This function also sets the boundary information in the
244 /// Node itself
245 //=========================================================
246 void Mesh::add_boundary_node(const unsigned &b, Node* const &node_pt)
247 {
248  //Tell the node that it's on boundary b.
249  //At this point, if the node is not a BoundaryNode, the function
250  //should throw an exception.
251  node_pt->add_to_boundary(b);
252 
253  //Get the size of the Boundary_node_pt vector
254  unsigned nbound_node=Boundary_node_pt[b].size();
255  bool node_already_on_this_boundary=false;
256  //Loop over the vector
257  for (unsigned n=0; n<nbound_node; n++)
258  {
259  // is the current node here already?
260  if (node_pt==Boundary_node_pt[b][n])
261  {
262  node_already_on_this_boundary=true;
263  }
264  }
265 
266  //Add the base node pointer to the vector if it's not there already
267  if (!node_already_on_this_boundary)
268  {
269  Boundary_node_pt[b].push_back(node_pt);
270  }
271 
272 }
273 
274 //=======================================================
275 /// Update nodal positions in response to changes in the domain shape.
276 /// Uses the FiniteElement::get_x(...) function for FiniteElements
277 /// and doesn't do anything for other element types.
278 /// If a MacroElement pointer has been set for a FiniteElement,
279 /// the MacroElement representation is used to update the
280 /// nodal positions; if not get_x(...) uses the FE interpolation
281 /// and thus leaves the nodal positions unchanged.
282 /// Virtual, so it can be overloaded by specific meshes,
283 /// such as AlgebraicMeshes or SpineMeshes.
284 /// Generally, this function updates the position of all nodes
285 /// in response to changes in the boundary position. For
286 /// SolidNodes it only applies the update to those SolidNodes
287 /// whose position is determined by the boundary position, unless
288 /// the bool flag is set to true.
289 //========================================================
290 void Mesh::node_update(const bool& update_all_solid_nodes)
291 {
292 #ifdef PARANOID
293 #ifdef OOMPH_HAS_MPI
294  //Paranoid check to throw an error if node update is called for elements
295  //with nonuniformly spaced nodes for which some masters are 'external'
296  for(unsigned long n=0;n<nnode();n++)
297  {
298  Node* nod_pt = Node_pt[n];
299  if (nod_pt->is_hanging())
300  {
301  //Loop over master nodes
302  unsigned nmaster=nod_pt->hanging_pt()->nmaster();
303  for (unsigned imaster=0;imaster<nmaster;imaster++)
304  {
305  //Get pointer to master node
306  Node* master_nod_pt = nod_pt->hanging_pt()->master_node_pt(imaster);
307 
308  //Get vector of all external halo nodes
310  get_external_halo_node_pt(external_halo_node_pt);
311 
312  //Search the external halo storage for this node
314  = std::find(external_halo_node_pt.begin(),
315  external_halo_node_pt.end(),
316  master_nod_pt);
317 
318  //Check if the node was found
319  if(it != external_halo_node_pt.end())
320  {
321  //Throw error becase node update won't work
322  //It's ok to throw an error here because this function is
323  //overloaded for Algebraic and MacroElementNodeUpdate
324  //Meshes. This is only a problem for meshes of ordinary
325  //nodes.
326  std::ostringstream err_stream;
327 
328  err_stream << "Calling node_update() for a mesh which contains"
329  << std::endl
330  << "master nodes which live in the external storage."
331  << std::endl
332  << "These nodes do not belong to elements on this"
333  << std::endl
334  << "processor and therefore cannot be updated locally."
335  << std::endl;
336 
337  throw OomphLibError(err_stream.str(),
338  OOMPH_CURRENT_FUNCTION,
339  OOMPH_EXCEPTION_LOCATION);
340  }
341  }
342  }
343  }
344  //If we get to here then none of the masters of any of the nodes in the
345  //mesh live in the external storage, so we'll be fine if we carry on.
346 #endif
347 #endif
348 
349  /// Local and global (Eulerian) coordinate
351  Vector<double> r;
352 
353  // NB: This repeats nodes a lot - surely it would be
354  // quicker to modify it so that it only does each node once,
355  // particularly in the update_all_solid_nodes=true case?
356 
357  // Loop over all elements
358  unsigned nel=nelement();
359  for (unsigned e=0;e<nel;e++)
360  {
361  // Try to cast to FiniteElement
362  FiniteElement* el_pt = dynamic_cast<FiniteElement*>(element_pt(e));
363 
364  // If it's a finite element we can proceed: FiniteElements have
365  // nodes and a get_x() function
366  if (el_pt!=0)
367  {
368  // Find out dimension of element = number of local coordinates
369  unsigned ndim_el=el_pt->dim();
370  s.resize(ndim_el);
371 
372  //Loop over nodal points
373  unsigned n_node=el_pt->nnode();
374  for (unsigned j=0;j<n_node;j++)
375  {
376 
377  // Get pointer to node
378  Node* nod_pt=el_pt->node_pt(j);
379 
380  // Get spatial dimension of node
381  unsigned ndim_node=nod_pt->ndim();
382  r.resize(ndim_node);
383 
384  // For non-hanging nodes
385  if (!(nod_pt->is_hanging()))
386  {
387  //Get the position of the node
388  el_pt->local_coordinate_of_node(j,s);
389 
390  // Get new position
391  el_pt->get_x(s,r);
392 
393  // Try to cast to SolidNode
394  SolidNode* solid_node_pt=dynamic_cast<SolidNode*>(nod_pt);
395 
396  // Loop over coordinate directions
397  for (unsigned i=0;i<ndim_node;i++)
398  {
399  // It's a SolidNode:
400  if (solid_node_pt!=0)
401  {
402  // only do it if explicitly requested!
403  if (update_all_solid_nodes)
404  {
405  solid_node_pt->x(i) = r[i];
406  }
407  }
408  // Not a SolidNode: Definitely update
409  else
410  {
411  nod_pt->x(i) = r[i];
412  }
413  }
414  }
415  }
416  }
417  }
418 
419  // Now update the external halo nodes before we adjust the positions of the
420  // hanging nodes incase any are masters of local nodes
421 #ifdef OOMPH_HAS_MPI
422  // Loop over all external halo nodes with other processors
423  // and update them
424  for (std::map<unsigned, Vector<Node*> >::iterator it=
425  External_halo_node_pt.begin();it!=External_halo_node_pt.end();it++)
426  {
427  // Get vector of external halo nodes
428  Vector<Node*> ext_halo_node_pt=(*it).second;
429  unsigned nnod=ext_halo_node_pt.size();
430  for (unsigned j=0;j<nnod;j++)
431  {
432  ext_halo_node_pt[j]->node_update();
433  }
434  }
435 #endif
436 
437  // Now loop over hanging nodes and adjust their position
438  // in line with their hanging node constraints
439  unsigned long n_node = nnode();
440  for(unsigned long n=0;n<n_node;n++)
441  {
442  Node* nod_pt = Node_pt[n];
443  if (nod_pt->is_hanging())
444  {
445  // Get spatial dimension of node
446  unsigned ndim_node=nod_pt->ndim();
447 
448  // Initialise
449  for (unsigned i=0;i<ndim_node;i++)
450  {
451  nod_pt->x(i)=0.0;
452  }
453 
454  //Loop over master nodes
455  unsigned nmaster=nod_pt->hanging_pt()->nmaster();
456  for (unsigned imaster=0;imaster<nmaster;imaster++)
457  {
458  // Loop over directions
459  for (unsigned i=0;i<ndim_node;i++)
460  {
461  nod_pt->x(i)+=nod_pt->hanging_pt()->
462  master_node_pt(imaster)->x(i)*
463  nod_pt->hanging_pt()->master_weight(imaster);
464  }
465  }
466  }
467  }
468 
469  // Loop over all nodes again and execute auxiliary node update
470  // function
471  for(unsigned long n=0;n<n_node;n++)
472  {
473  Node_pt[n]->perform_auxiliary_node_update_fct();
474  }
475 
476 }
477 
478 
479 //=======================================================
480 /// Reorder nodes in the order in which they are
481 /// encountered when stepping through the elements
482 //========================================================
483  void Mesh::reorder_nodes(const bool& use_old_ordering)
484  {
485  Vector<Node*> reordering;
486  get_node_reordering(reordering, use_old_ordering);
487 
488  unsigned n_node = nnode();
489  for(unsigned i=0; i<n_node; i++)
490  {
491  node_pt(i) = reordering[i];
492  }
493  }
494 
495 //=======================================================
496 /// Get a vector of the nodes in the order in which they are encountered
497 /// when stepping through the elements (similar to reorder_nodes() but
498 /// without changing the mesh's node vector).
499 //========================================================
501  const bool& use_old_ordering) const
502  {
503  if(use_old_ordering)
504  {
505  // Setup map to check if nodes have been done yet
506  std::map<Node*,bool> done;
507 
508  // Loop over all nodes
509  unsigned nnod=nnode();
510 
511  // Initialise the vector
512  reordering.assign(nnod,0);
513 
514  // Return immediately if there are no nodes: Note assumption:
515  // Either all the elements' nodes stored here or none.
516  // If only a subset is stored in the Node_pt vector we'll get
517  // a range checking error below (only if run with range, checking,
518  // of course).
519  if (nnod==0) return;
520  for (unsigned j=0;j<nnod;j++)
521  {
522  done[node_pt(j)]=false;
523  }
524 
525  // Initialise counter for number of nodes
526  unsigned long count=0;
527 
528  // Loop over all elements
529  unsigned nel=nelement();
530  for (unsigned e=0;e<nel;e++)
531  {
532  // Loop over nodes in element
534  unsigned nnod=el_pt->nnode();
535  for (unsigned j=0;j<nnod;j++)
536  {
537  Node* nod_pt=el_pt->node_pt(j);
538  // Has node been done yet?
539  if (!done[nod_pt])
540  {
541  // Insert into node vector. NOTE: If you get a seg fault/range checking
542  // error here then you probably haven't added all the elements' nodes
543  // to the Node_pt vector -- this is most likely to arise in the
544  // case of meshes of face elements (though they usually don't
545  // store the nodes at all so if you have any problems here there's
546  // something unusual/not quite right in any case... For this
547  // reason we don't range check here by default (not even under
548  // paranoia) but force you turn on proper (costly) range checking
549  // to track this down...
550  reordering[count]=nod_pt;
551  done[nod_pt]=true;
552  // Increase counter
553  count++;
554  }
555  }
556  }
557 
558  // Sanity check
559  if (count!=nnod)
560  {
561  throw OomphLibError(
562  "Trouble: Number of nodes hasn't stayed constant during reordering!\n",
563  OOMPH_CURRENT_FUNCTION,
564  OOMPH_EXCEPTION_LOCATION);
565  }
566 
567  }
568  else
569  {
570  // Copy node vector out
571  unsigned n_node = nnode();
572  reordering.resize(n_node);
573  for(unsigned i=0; i<n_node; i++)
574  {
575  reordering[i] = node_pt(i);
576  }
577 
578  // and sort it
579  std::sort(reordering.begin(), reordering.end(),
581  }
582  }
583 
584 
585 //========================================================
586 /// Virtual Destructor to clean up all memory
587 //========================================================
589 {
590  //Free the nodes first
591  //Loop over the nodes in reverse order
592  unsigned long Node_pt_range = Node_pt.size();
593  for(unsigned long i=Node_pt_range;i>0;i--)
594  {
595  delete Node_pt[i-1]; Node_pt[i-1] = 0;
596  }
597 
598  //Free the elements
599  //Loop over the elements in reverse order
600  unsigned long Element_pt_range = Element_pt.size();
601  for(unsigned long i=Element_pt_range;i>0;i--)
602  {
603  delete Element_pt[i-1]; Element_pt[i-1] = 0;
604  }
605 
606  // Wipe the storage for all externally-based elements and delete halos
608 
609 }
610 
611 //========================================================
612 /// Assign (global) equation numbers to the nodes
613 //========================================================
615 {
616  //Find out the current number of equations
617  unsigned long equation_number=Dof_pt.size();
618 
619  //Loop over the nodes and call their assigment functions
620  unsigned long nnod = Node_pt.size();
621 
622  for(unsigned long i=0;i<nnod;i++)
623  {
624  Node_pt[i]->assign_eqn_numbers(equation_number,Dof_pt);
625  }
626 
627  //Loop over the elements and number their internals
628  unsigned long nel = Element_pt.size();
629  for(unsigned long i=0;i<nel;i++)
630  {
631  Element_pt[i]->assign_internal_eqn_numbers(equation_number,Dof_pt);
632  }
633 
634  //Return the total number of equations
635  return(equation_number);
636 }
637 
638 //========================================================
639 /// \short Function to describe the dofs of the Mesh. The ostream
640 /// specifies the output stream to which the description
641 /// is written; the string stores the currently
642 /// assembled output that is ultimately written to the
643 /// output stream by Data::describe_dofs(...); it is typically
644 /// built up incrementally as we descend through the
645 /// call hierarchy of this function when called from
646 /// Problem::describe_dofs(...)
647 //========================================================
648 void Mesh::describe_dofs(std::ostream& out,
649  const std::string& current_string) const
650 {
651  //Loop over the nodes and call their classification functions
652  unsigned long nnod = Node_pt.size();
653  for(unsigned long i=0;i<nnod;i++)
654  {
655  std::stringstream conversion;
656  conversion << " of Node " << i << current_string;
657  std::string in(conversion.str());
658  Node_pt[i]->describe_dofs(out,in);
659  }
660 
661  //Loop over the elements and classify.
662  unsigned long nel = Element_pt.size();
663  for(unsigned long i=0;i<nel;i++)
664  {
665  std::stringstream conversion;
666  conversion <<" in Element "<<i<<" ["<<typeid(*Element_pt[i]).name()<<"] "
667  <<current_string;
668  std::string in(conversion.str());
669  Element_pt[i]->describe_dofs(out,in);
670  }
671 }
672 
673 //========================================================
674 /// \short Function to describe the local dofs of the elements. The ostream
675 /// specifies the output stream to which the description
676 /// is written; the string stores the currently
677 /// assembled output that is ultimately written to the
678 /// output stream by Data::describe_dofs(...); it is typically
679 /// built up incrementally as we descend through the
680 /// call hierarchy of this function when called from
681 /// Problem::describe_dofs(...)
682 //========================================================
683 void Mesh::describe_local_dofs(std::ostream& out,
684  const std::string& current_string) const
685 {
686  //Now loop over the elements and describe local dofs
687  unsigned long nel = Element_pt.size();
688  for(unsigned long i=0;i<nel;i++)
689  {
690  std::stringstream conversion;
691  conversion <<" in Element"<<i<<" ["<<typeid(*Element_pt[i]).name()<<"] "
692  << current_string;
693  std::string in(conversion.str());
694  Element_pt[i]->describe_local_dofs(out,in);
695  }
696 }
697 
698 
699 //========================================================
700 /// Assign local equation numbers in all elements
701 //========================================================
702 void Mesh::assign_local_eqn_numbers(const bool &store_local_dof_pt)
703 {
704  //Now loop over the elements and assign local equation numbers
705  unsigned long Element_pt_range = Element_pt.size();
706  for(unsigned long i=0;i<Element_pt_range;i++)
707  {
708  Element_pt[i]->assign_local_eqn_numbers(store_local_dof_pt);
709  }
710 }
711 
712 //========================================================
713 /// Self-test: Check elements and nodes. Return 0 for OK
714 //========================================================
715 unsigned Mesh::self_test()
716 {
717 
718  // Initialise
719  bool passed=true;
720 
721  // Check the mesh for repeated nodes (issues its own error message)
722  if (0!=check_for_repeated_nodes()) passed=false;
723 
724 // hierher -- re-enable once problem with Hermite elements has been resolved.
725 // // Check if there are any inverted elements
726 // bool mesh_has_inverted_elements=false;
727 // check_inverted_elements(mesh_has_inverted_elements);
728 // if (mesh_has_inverted_elements)
729 // {
730 // passed=false;
731 // oomph_info << "\n ERROR: Mesh has inverted elements\n"
732 // << " Run Mesh::check_inverted_elements(...) with"
733 // << " with output stream to find out which elements are"
734 // << " inverted.\n";
735 // }
736 
737  //Loop over the elements, check for duplicates and do self test
738  std::set<GeneralisedElement*> element_set_pt;
739  unsigned long Element_pt_range = Element_pt.size();
740  for(unsigned long i=0;i<Element_pt_range;i++)
741  {
742  if (Element_pt[i]->self_test()!=0)
743  {
744  passed=false;
745  oomph_info << "\n ERROR: Failed Element::self_test() for element i="
746  << i << std::endl;
747  }
748  // Add to set (which ignores duplicates):
749  element_set_pt.insert(Element_pt[i]);
750  }
751 
752  // Check for duplicates:
753  if (element_set_pt.size()!=Element_pt_range)
754  {
755  oomph_info << "ERROR: " << Element_pt_range-element_set_pt.size()
756  << " duplicate elements were encountered in mesh!" << std::endl;
757  passed=false;
758  }
759 
760 
761  //Loop over the nodes, check for duplicates and do self test
762  std::set<Node*> node_set_pt;
763  unsigned long Node_pt_range = Node_pt.size();
764  for(unsigned long i=0;i<Node_pt_range;i++)
765  {
766  if (Node_pt[i]->self_test()!=0)
767  {
768  passed=false;
769  oomph_info << "\n ERROR: Failed Node::self_test() for node i="
770  << i << std::endl;
771  }
772  // Add to set (which ignores duplicates):
773  node_set_pt.insert(Node_pt[i]);
774  }
775 
776  // Check for duplicates:
777  if (node_set_pt.size()!=Node_pt_range)
778  {
779  oomph_info << "ERROR: " << Node_pt_range-node_set_pt.size()
780  << " duplicate nodes were encountered in mesh!" << std::endl;
781  passed=false;
782  }
783 
784  // Return verdict
785  if (passed) {return 0;}
786  else {return 1;}
787 
788 }
789 
790 
791 
792 //========================================================
793 /// Check for inverted elements and report outcome
794  /// in boolean variable. This visits all elements at their
795  /// integration points and checks if the Jacobian of the
796  /// mapping between local and global coordinates is positive --
797  /// using the same test that would be carried out (but only in PARANOID
798  /// mode) during the assembly of the elements' Jacobian matrices.
799  /// Inverted elements are output in inverted_element_file (if the
800  /// stream is open).
801 //========================================================
802  void Mesh::check_inverted_elements(bool& mesh_has_inverted_elements,
803  std::ofstream& inverted_element_file)
804  {
805  // Initialise flag
806  mesh_has_inverted_elements=false;
807 
808  // Suppress output while checking for inverted elements
809  bool backup=
812 
813  // Loop over all elements
814  unsigned nelem=nelement();
815  for (unsigned e=0;e<nelem;e++)
816  {
818 
819  // Only check for finite elements
820  if (el_pt!=0)
821  {
822 
823  //Find out number of nodes and local coordinates in the element
824  unsigned n_node = el_pt->nnode();
825  unsigned n_dim=el_pt->dim();
826  unsigned ndim_node=el_pt->nodal_dimension();
827 
828  // Can't check Jacobian for elements in which nodal and elementa
829  // dimensions don't match
830  if (n_dim==ndim_node)
831  {
832 
833  //Set up memory for the shape and test function and local coord
834  Shape psi(n_node);
835  DShape dpsidx(n_node,n_dim);
836  Vector<double> s(n_dim);
837 
838  // Initialise element-level test
839  bool is_inverted=false;
840 
841  unsigned n_intpt=el_pt->integral_pt()->nweight();
842  for(unsigned ipt=0;ipt<n_intpt;ipt++)
843  {
844  // Local coordinates
845  for(unsigned i=0;i<n_dim;i++)
846  {
847  s[i] = el_pt->integral_pt()->knot(ipt,i);
848  }
849 
850  double J=0;
851  // Dummy assignment to keep gcc from complaining about
852  // "set but unused".
853  J+=0.0;
854  try
855  {
856  //Call the derivatives of the shape functions and Jacobian
857  J=el_pt->dshape_eulerian(s,psi,dpsidx);
858 
859  // If code is compiled without PARANOID setting, the
860  // above call will simply return the negative Jacobian
861  // without failing, so we need to check manually
862 #ifndef PARANOID
863  try
864  {
865  el_pt->check_jacobian(J);
866  }
867  catch(OomphLibQuietException& error)
868  {
869  is_inverted=true;
870  }
871 #endif
872  }
873  catch(OomphLibQuietException& error)
874  {
875  is_inverted=true;
876  }
877  }
878  if (is_inverted)
879  {
880  mesh_has_inverted_elements=true;
881  if (inverted_element_file.is_open())
882  {
883  el_pt->output(inverted_element_file);
884  }
885  }
886  }
887  }
888  }
889  // Reset
891  }
892 
893 
894 
895 //========================================================
896 /// Nodes that have been marked as obsolete are removed
897 /// from the mesh and the its boundaries. Returns vector
898 /// of pointers to deleted nodes.
899 //========================================================
901 {
902 
903  // Only copy the 'live' nodes across to new mesh
904  //----------------------------------------------
905 
906  // New Vector of pointers to nodes
907  Vector<Node *> new_node_pt;
908  Vector<Node *> deleted_node_pt;
909 
910  // Loop over all nodes in mesh
911  unsigned long n_node = nnode();
912  for(unsigned long n=0;n<n_node;n++)
913  {
914  // If the node still exists: Copy across
915  if(!(Node_pt[n]->is_obsolete()))
916  {
917  new_node_pt.push_back(Node_pt[n]);
918  }
919  // Otherwise the Node is gone:
920  // Delete it for good if it does not lie on a boundary
921  // (if it lives on a boundary we have to remove it from
922  // the boundary lookup schemes below)
923  else
924  {
925  if(!(Node_pt[n]->is_on_boundary()))
926  {
927  deleted_node_pt.push_back(Node_pt[n]);
928  delete Node_pt[n];
929  Node_pt[n]=0;
930  }
931  }
932  }
933 
934  // Now update old vector by setting it equal to the new vector
935  Node_pt = new_node_pt;
936 
937 
938  // Boundaries
939  //-----------
940 
941  // Only copy the 'live' nodes into new boundary node arrays
942  //---------------------------------------------------------
943  //Loop over the boundaries
944  unsigned num_bound = nboundary();
945  for(unsigned ibound=0;ibound<num_bound;ibound++)
946  {
947  // New Vector of pointers to existent boundary nodes
948  Vector<Node *> new_boundary_node_pt;
949 
950  //Loop over the boundary nodes
951  unsigned long Nboundary_node = Boundary_node_pt[ibound].size();
952 
953  // Reserve contiguous memory for new vector of pointers
954  // Must be equal in size to the number of nodes or less
955  new_boundary_node_pt.reserve(Nboundary_node);
956 
957  for(unsigned long n=0;n<Nboundary_node;n++)
958  {
959  // If node still exists: Copy across
960  if (!(Boundary_node_pt[ibound][n]->is_obsolete()))
961  {
962  new_boundary_node_pt.push_back(Boundary_node_pt[ibound][n]);
963  }
964  // Otherwise Node is gone: Delete it for good
965  else
966  {
967  //The node may lie on multiple boundaries, so remove the node
968  //from the current boundary
969  Boundary_node_pt[ibound][n]->remove_from_boundary(ibound);
970 
971  //Now if the node is no longer on any boundaries, delete it
972  if(!Boundary_node_pt[ibound][n]->is_on_boundary())
973  {
974  deleted_node_pt.push_back(
975  dynamic_cast<Node*>(Boundary_node_pt[ibound][n]));
976 
977  delete Boundary_node_pt[ibound][n];
978  }
979  }
980  }
981 
982  //Update the old vector by setting it equal to the new vector
983  Boundary_node_pt[ibound] = new_boundary_node_pt;
984 
985  } //End of loop over boundaries
986 
987  // Tell us who you deleted
988  return deleted_node_pt;
989 }
990 
991 
992 
993 
994 //========================================================
995 /// Output function for the mesh boundaries
996 ///
997 /// Loop over all boundaries and dump out the coordinates
998 /// of the points on the boundary (in individual tecplot
999 /// zones)
1000 //========================================================
1001 void Mesh::output_boundaries(std::ostream &outfile)
1002 {
1003  //Loop over the boundaries
1004  unsigned num_bound = nboundary();
1005  for(unsigned long ibound=0;ibound<num_bound;ibound++)
1006  {
1007  unsigned nnod=Boundary_node_pt[ibound].size();
1008  if (nnod>0)
1009  {
1010  outfile << "ZONE T=\"boundary" << ibound << "\"\n";
1011 
1012  for (unsigned inod=0;inod<nnod;inod++)
1013  {
1014  Boundary_node_pt[ibound][inod]->output(outfile);
1015  }
1016  }
1017  }
1018 }
1019 
1020 //========================================================
1021 /// Output the nodes on the boundaries and their
1022 /// respective boundary coordinates(into separate tecplot
1023 /// zones)
1024 //========================================================
1025 void Mesh::output_boundaries_coordinates(std::ostream &outfile)
1026 {
1027  //Loop over the boundaries
1028  const unsigned nbound = nboundary();
1029  for(unsigned ibound = 0; ibound < nbound; ibound++)
1030  {
1031  unsigned nelement = Boundary_element_pt[ibound].size();
1032  if (nelement>0)
1033  {
1034  output_boundary_coordinates(ibound, outfile);
1035  } // if (nelement > 0)
1036  } // for (ibound < nbound)
1037 }
1038 
1039 
1040 //===================================================================
1041 /// Dump function for the mesh class.
1042 /// Loop over all nodes and elements and dump them
1043 //===================================================================
1044  void Mesh::dump(std::ofstream &dump_file, const bool& use_old_ordering) const
1045 {
1046 
1047  // Get a reordering of the nodes so that the dump file is in a standard
1048  // ordering regardless of the sequence of mesh refinements etc.
1049  Vector<Node*> reordering;
1050  this->get_node_reordering(reordering, use_old_ordering);
1051 
1052  // Find number of nodes
1053  unsigned long Node_pt_range = this->nnode();
1054 
1055  // Doc # of nodes
1056  dump_file << Node_pt_range << " # number of nodes " << std::endl;
1057 
1058  //Loop over all the nodes and dump their data
1059  for(unsigned nd=0; nd<Node_pt_range; nd++)
1060  {
1061  reordering[nd]->dump(dump_file);
1062  }
1063 
1064  // Loop over elements and deal with internal data
1065  unsigned n_element = this->nelement();
1066  for(unsigned e=0;e<n_element;e++)
1067  {
1068  GeneralisedElement* el_pt = this->element_pt(e);
1069  unsigned n_internal = el_pt->ninternal_data();
1070  if(n_internal > 0)
1071  {
1072  dump_file << n_internal
1073  << " # number of internal Data items in element "
1074  << e << std::endl;
1075  for(unsigned i=0;i<n_internal;i++)
1076  {
1077  el_pt->internal_data_pt(i)->dump(dump_file);
1078  }
1079  }
1080  }
1081 }
1082 
1083 
1084 //=======================================================
1085 /// Read solution from restart file
1086 //=======================================================
1087 void Mesh::read(std::ifstream &restart_file)
1088 {
1089  std::string input_string;
1090 
1091  // Reorder the nodes within the mesh's node vector
1092  // to establish a standard ordering regardless of the sequence
1093  // of mesh refinements etc
1094  this->reorder_nodes();
1095 
1096  //Read nodes
1097 
1098  // Find number of nodes
1099  unsigned long n_node = this->nnode();
1100 
1101  // Read line up to termination sign
1102  getline(restart_file,input_string,'#');
1103 
1104  // Ignore rest of line
1105  restart_file.ignore(80,'\n');
1106 
1107  // Check # of nodes:
1108  unsigned long check_n_node=atoi(input_string.c_str());
1109  if (check_n_node!=n_node)
1110  {
1111  std::ostringstream error_stream;
1112  error_stream << "The number of nodes allocated " << n_node
1113  << " is not the same as specified in the restart file "
1114  << check_n_node << std::endl;
1115 
1116  throw OomphLibError(error_stream.str(),
1117  OOMPH_CURRENT_FUNCTION,
1118  OOMPH_EXCEPTION_LOCATION);
1119  }
1120 
1121  //Loop over the nodes
1122  for(unsigned long n=0;n<n_node;n++)
1123  {
1124  /// Try to cast to elastic node
1125  SolidNode* el_node_pt=dynamic_cast<SolidNode*>(
1126  this->node_pt(n));
1127  if (el_node_pt!=0)
1128  {
1129  el_node_pt->read(restart_file);
1130  }
1131  else
1132  {
1133  this->node_pt(n)->read(restart_file);
1134  }
1135  }
1136 
1137  // Read internal data of elements:
1138  //--------------------------------
1139  // Loop over elements and deal with internal data
1140  unsigned n_element = this->nelement();
1141  for (unsigned e=0;e<n_element;e++)
1142  {
1143  GeneralisedElement* el_pt = this->element_pt(e);
1144  unsigned n_internal=el_pt->ninternal_data();
1145  if (n_internal>0)
1146  {
1147  // Read line up to termination sign
1148  getline(restart_file,input_string,'#');
1149 
1150  // Ignore rest of line
1151  restart_file.ignore(80,'\n');
1152 
1153  // Check # of internals :
1154  unsigned long check_n_internal=atoi(input_string.c_str());
1155  if (check_n_internal!=n_internal)
1156  {
1157  std::ostringstream error_stream;
1158  error_stream << "The number of internal data " << n_internal
1159  << " is not the same as specified in the restart file "
1160  << check_n_internal << std::endl;
1161 
1162  throw OomphLibError(error_stream.str(),
1163  OOMPH_CURRENT_FUNCTION,
1164  OOMPH_EXCEPTION_LOCATION);
1165  }
1166 
1167  for (unsigned i=0;i<n_internal;i++)
1168  {
1169  el_pt->internal_data_pt(i)->read(restart_file);
1170  }
1171  }
1172  }
1173 }
1174 
1175 
1176 //========================================================
1177 /// Output in paraview format into specified file.
1178 ///
1179 /// Breaks up each element into sub-elements for plotting
1180 /// purposes. We assume that all elements are of the same
1181 /// type (fct will break (in paranoid mode) if paraview
1182 /// output fcts of the elements are inconsistent).
1183 //========================================================
1184 void Mesh::output_paraview(std::ofstream &file_out,
1185  const unsigned &nplot) const
1186 {
1187 
1188  // Change the scientific format so that E is used rather than e
1189  file_out.setf(std::ios_base::uppercase);
1190 
1191  // Decide how many elements there are to be plotted
1192  unsigned long number_of_elements=this->Element_pt.size();
1193 
1194  // Cast to finite element and return if cast fails.
1195  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(0));
1196 
1197 #ifdef PARANOID
1198  if (fe_pt==0)
1199  {
1200  throw OomphLibError(
1201  "Recast for FiniteElement failed for element 0!\n",
1202  OOMPH_CURRENT_FUNCTION,
1203  OOMPH_EXCEPTION_LOCATION);
1204  }
1205 #endif
1206 
1207 
1208 #ifdef PARANOID
1209  // Check if all elements have the same number of degrees of freedom,
1210  // if they don't, paraview will break
1211  unsigned el_zero_ndof=fe_pt->nscalar_paraview();
1212  for(unsigned i=1;i<number_of_elements;i++)
1213  {
1214  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1215  unsigned el_i_ndof=fe_pt->nscalar_paraview();
1216  if(el_zero_ndof!=el_i_ndof)
1217  {
1218  std::stringstream error_stream;
1219  error_stream
1220  << "Element " << i << " has different number of degrees of freedom\n"
1221  << "than from previous elements, Paraview cannot handle this.\n"
1222  << "We suggest that the problem is broken up into submeshes instead."
1223  << std::endl;
1224  throw OomphLibError(
1225  error_stream.str(),
1226  OOMPH_CURRENT_FUNCTION,
1227  OOMPH_EXCEPTION_LOCATION);
1228  }
1229  }
1230 #endif
1231 
1232  // Make variables to hold the number of nodes and elements
1233  unsigned long number_of_nodes=0;
1234  unsigned long total_number_of_elements=0;
1235 
1236  // Loop over all the elements to find total number of plot points
1237  for(unsigned i=0;i<number_of_elements;i++)
1238  {
1239  // Cast to FiniteElement and (in paranoid mode) check
1240  // if cast has failed.
1241  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1242 
1243 #ifdef PARANOID
1244  if (fe_pt==0)
1245  {
1246  std::stringstream error_stream;
1247  error_stream
1248  << "Recast for element " << i << " failed" << std::endl;
1249  throw OomphLibError(
1250  error_stream.str(),
1251  OOMPH_CURRENT_FUNCTION,
1252  OOMPH_EXCEPTION_LOCATION);
1253  }
1254 #endif
1255 
1256  number_of_nodes+=fe_pt->nplot_points_paraview(nplot);
1257  total_number_of_elements+=fe_pt->nsub_elements_paraview(nplot);
1258 
1259  }
1260 
1261 
1262  // File Declaration
1263  //------------------
1264 
1265  // Insert the necessary lines plus header of file, and
1266  // number of nodes and elements
1267  file_out
1268  << "<?xml version=\"1.0\"?>\n"
1269  << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
1270  << "byte_order=\"LittleEndian\">\n"
1271  << "<UnstructuredGrid>\n"
1272  << "<Piece NumberOfPoints=\""
1273  << number_of_nodes
1274  << "\" NumberOfCells=\""
1275  << total_number_of_elements
1276  <<"\">\n";
1277 
1278 
1279  // Point Data
1280  //-----------
1281 
1282  // Check the number of degrees of freedom
1283  unsigned ndof = fe_pt->nscalar_paraview();
1284 
1285  // Point data is going in here
1286  file_out << "<PointData ";
1287 
1288  // Insert just the first scalar name, since paraview reads everything
1289  // else after that as being of the same type. Get information from
1290  // first element.
1291  file_out << "Scalars=\""
1292  << fe_pt->scalar_name_paraview(0)
1293  << "\">\n";
1294 
1295  // Loop over i scalar fields and j number of elements
1296  for(unsigned i=0;i<ndof;i++)
1297  {
1298  file_out << "<DataArray type=\"Float32\" "
1299  << "Name=\""
1300  << fe_pt->scalar_name_paraview(i)
1301  << "\" "
1302  << "format=\"ascii\""
1303  << ">\n";
1304 
1305  for(unsigned j=0;j<number_of_elements;j++)
1306  {
1307  // Cast to FiniteElement and (in paranoid mode) check
1308  // if cast has failed.
1309  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(j));
1310 
1311 #ifdef PARANOID
1312  if (fe_pt==0)
1313  {
1314  std::stringstream error_stream;
1315  error_stream
1316  << "Recast for element " << j << " failed" << std::endl;
1317  throw OomphLibError(
1318  error_stream.str(),
1319  OOMPH_CURRENT_FUNCTION,
1320  OOMPH_EXCEPTION_LOCATION);
1321  }
1322 #endif
1323 
1324  fe_pt->scalar_value_paraview(file_out,i,nplot);
1325  }
1326 
1327  // Close of the DataArray
1328  file_out << "</DataArray>\n";
1329  }
1330 
1331  // Close off the PointData set
1332  file_out << "</PointData>\n";
1333 
1334 
1335  // Geometric Points
1336  //------------------
1337 
1338  file_out
1339  << "<Points>\n"
1340  << "<DataArray type=\"Float32\""
1341  << " NumberOfComponents=\""
1342  // This always has to be 3 for an unstructured grid
1343  << 3 << "\" "
1344  << "format=\"ascii\">\n";
1345 
1346  // Loop over all the elements to print their plot points
1347  for(unsigned i=0;i<number_of_elements;i++)
1348  {
1349  // Cast to FiniteElement and (in paranoid mode) check
1350  // if cast has failed.
1351  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1352 
1353 #ifdef PARANOID
1354  if (fe_pt==0)
1355  {
1356  std::stringstream error_stream;
1357  error_stream
1358  << "Recast for element " << i << " faild" << std::endl;
1359  throw OomphLibError(
1360  error_stream.str(),
1361  OOMPH_CURRENT_FUNCTION,
1362  OOMPH_EXCEPTION_LOCATION);
1363  }
1364 #endif
1365 
1366  fe_pt->output_paraview(file_out,nplot);
1367  }
1368 
1369  file_out
1370  << "</DataArray>\n"
1371  << "</Points>\n";
1372 
1373 
1374  // Cells
1375  //-------
1376 
1377  file_out
1378  << "<Cells>\n"
1379  << "<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";
1380 
1381  // Make counter for keeping track of all the local elements,
1382  // because Paraview requires global coordinates
1383  unsigned counter=0;
1384 
1385  // Write connectivity with the local elements
1386  for(unsigned i=0;i<number_of_elements;i++)
1387  {
1388  // Cast to FiniteElement and (in paranoid mode) check
1389  // if cast has failed.
1390  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1391 
1392 #ifdef PARANOID
1393  if (fe_pt==0)
1394  {
1395  std::stringstream error_stream;
1396  error_stream
1397  << "Recast for element " << i << " faild" << std::endl;
1398  throw OomphLibError(
1399  error_stream.str(),
1400  OOMPH_CURRENT_FUNCTION,
1401  OOMPH_EXCEPTION_LOCATION);
1402  }
1403 #endif
1404  fe_pt->write_paraview_output_offset_information(file_out,nplot,counter);
1405  }
1406 
1407  file_out << "</DataArray>\n"
1408  << "<DataArray type=\"Int32\" "
1409  << "Name=\"offsets\" format=\"ascii\">\n";
1410 
1411  // Make variable that holds the current offset number
1412  unsigned offset_sum=0;
1413 
1414  // Write the offset for the specific elements
1415  for(unsigned i=0;i<number_of_elements;i++)
1416  {
1417  // Cast to FiniteElement and (in paranoid mode) check
1418  // if cast has failed.
1419  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1420 
1421 #ifdef PARANOID
1422  if (fe_pt==0)
1423  {
1424  std::stringstream error_stream;
1425  error_stream
1426  << "Recast for element " << i << " failed" << std::endl;
1427  throw OomphLibError(
1428  error_stream.str(),
1429  OOMPH_CURRENT_FUNCTION,
1430  OOMPH_EXCEPTION_LOCATION);
1431  }
1432 #endif
1433  fe_pt->write_paraview_offsets(file_out,nplot,offset_sum);
1434  }
1435 
1436  file_out <<"</DataArray>\n"
1437  <<"<DataArray type=\"UInt8\" Name=\"types\">\n";
1438 
1439  // Loop over all elements to get the type that they have
1440  for(unsigned i=0;i<number_of_elements;i++)
1441  {
1442  // Cast to FiniteElement and (in paranoid mode) check
1443  // if cast has failed.
1444  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1445 
1446 #ifdef PARANOID
1447  if (fe_pt==0)
1448  {
1449  std::stringstream error_stream;
1450  error_stream
1451  << "Recast for element " << i << " failed" << std::endl;
1452  throw OomphLibError(
1453  error_stream.str(),
1454  OOMPH_CURRENT_FUNCTION,
1455  OOMPH_EXCEPTION_LOCATION);
1456  }
1457 #endif
1458 
1459  fe_pt->write_paraview_type(file_out,nplot);
1460  }
1461 
1462  file_out <<"</DataArray>\n"
1463  <<"</Cells>\n";
1464 
1465 
1466  // File Closure
1467  //-------------
1468  file_out <<"</Piece>\n"
1469  <<"</UnstructuredGrid>\n"
1470  <<"</VTKFile>";
1471 }
1472 
1473 
1474 //========================================================
1475 /// Output in paraview format into specified file.
1476 ///
1477 /// Breaks up each element into sub-elements for plotting
1478 /// purposes. We assume that all elements are of the same
1479 /// type (fct will break (in paranoid mode) if paraview
1480 /// output fcts of the elements are inconsistent).
1481 //========================================================
1482  void Mesh::output_fct_paraview(std::ofstream &file_out,
1483  const unsigned &nplot,
1484  FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
1485 {
1486 
1487  // Change the scientific format so that E is used rather than e
1488  file_out.setf(std::ios_base::uppercase);
1489 
1490  // Decide how many elements there are to be plotted
1491  unsigned long number_of_elements=this->Element_pt.size();
1492 
1493  // Cast to finite element and return if cast fails.
1494  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(0));
1495 
1496 #ifdef PARANOID
1497  if (fe_pt==0)
1498  {
1499  throw OomphLibError(
1500  "Recast for FiniteElement failed for element 0!\n",
1501  OOMPH_CURRENT_FUNCTION,
1502  OOMPH_EXCEPTION_LOCATION);
1503  }
1504 #endif
1505 
1506 
1507 #ifdef PARANOID
1508  // Check if all elements have the same number of degrees of freedom,
1509  // if they don't, paraview will break
1510  unsigned el_zero_ndof=fe_pt->nscalar_paraview();
1511  for(unsigned i=1;i<number_of_elements;i++)
1512  {
1513  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1514  unsigned el_i_ndof=fe_pt->nscalar_paraview();
1515  if(el_zero_ndof!=el_i_ndof)
1516  {
1517  std::stringstream error_stream;
1518  error_stream
1519  << "Element " << i << " has different number of degrees of freedom\n"
1520  << "than from previous elements, Paraview cannot handle this.\n"
1521  << "We suggest that the problem is broken up into submeshes instead."
1522  << std::endl;
1523  throw OomphLibError(
1524  error_stream.str(),
1525  OOMPH_CURRENT_FUNCTION,
1526  OOMPH_EXCEPTION_LOCATION);
1527  }
1528  }
1529 #endif
1530 
1531  // Make variables to hold the number of nodes and elements
1532  unsigned long number_of_nodes=0;
1533  unsigned long total_number_of_elements=0;
1534 
1535  // Loop over all the elements to find total number of plot points
1536  for(unsigned i=0;i<number_of_elements;i++)
1537  {
1538  // Cast to FiniteElement and (in paranoid mode) check
1539  // if cast has failed.
1540  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1541 
1542 #ifdef PARANOID
1543  if (fe_pt==0)
1544  {
1545  std::stringstream error_stream;
1546  error_stream
1547  << "Recast for element " << i << " failed" << std::endl;
1548  throw OomphLibError(
1549  error_stream.str(),
1550  OOMPH_CURRENT_FUNCTION,
1551  OOMPH_EXCEPTION_LOCATION);
1552  }
1553 #endif
1554 
1555  number_of_nodes+=fe_pt->nplot_points_paraview(nplot);
1556  total_number_of_elements+=fe_pt->nsub_elements_paraview(nplot);
1557 
1558  }
1559 
1560 
1561  // File Declaration
1562  //------------------
1563 
1564  // Insert the necessary lines plus header of file, and
1565  // number of nodes and elements
1566  file_out
1567  << "<?xml version=\"1.0\"?>\n"
1568  << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
1569  << "byte_order=\"LittleEndian\">\n"
1570  << "<UnstructuredGrid>\n"
1571  << "<Piece NumberOfPoints=\""
1572  << number_of_nodes
1573  << "\" NumberOfCells=\""
1574  << total_number_of_elements
1575  <<"\">\n";
1576 
1577 
1578  // Point Data
1579  //-----------
1580 
1581  // Check the number of degrees of freedom
1582  unsigned ndof = fe_pt->nscalar_paraview();
1583 
1584  // Point data is going in here
1585  file_out << "<PointData ";
1586 
1587  // Insert just the first scalar name, since paraview reads everything
1588  // else after that as being of the same type. Get information from
1589  // first element.
1590  file_out << "Scalars=\""
1591  << fe_pt->scalar_name_paraview(0)
1592  << "\">\n";
1593 
1594  // Loop over i scalar fields and j number of elements
1595  for(unsigned i=0;i<ndof;i++)
1596  {
1597  file_out << "<DataArray type=\"Float32\" "
1598  << "Name=\""
1599  << fe_pt->scalar_name_paraview(i)
1600  << "\" "
1601  << "format=\"ascii\""
1602  << ">\n";
1603 
1604  for(unsigned j=0;j<number_of_elements;j++)
1605  {
1606  // Cast to FiniteElement and (in paranoid mode) check
1607  // if cast has failed.
1608  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(j));
1609 
1610 #ifdef PARANOID
1611  if (fe_pt==0)
1612  {
1613  std::stringstream error_stream;
1614  error_stream
1615  << "Recast for element " << j << " failed" << std::endl;
1616  throw OomphLibError(
1617  error_stream.str(),
1618  OOMPH_CURRENT_FUNCTION,
1619  OOMPH_EXCEPTION_LOCATION);
1620  }
1621 #endif
1622 
1623  fe_pt->scalar_value_fct_paraview(file_out,i,nplot,exact_soln_pt);
1624  }
1625 
1626  // Close of the DataArray
1627  file_out << "</DataArray>\n";
1628  }
1629 
1630  // Close off the PointData set
1631  file_out << "</PointData>\n";
1632 
1633 
1634  // Geometric Points
1635  //------------------
1636 
1637  file_out
1638  << "<Points>\n"
1639  << "<DataArray type=\"Float32\""
1640  << " NumberOfComponents=\""
1641  // This always has to be 3 for an unstructured grid
1642  << 3 << "\" "
1643  << "format=\"ascii\">\n";
1644 
1645  // Loop over all the elements to print their plot points
1646  for(unsigned i=0;i<number_of_elements;i++)
1647  {
1648  // Cast to FiniteElement and (in paranoid mode) check
1649  // if cast has failed.
1650  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1651 
1652 #ifdef PARANOID
1653  if (fe_pt==0)
1654  {
1655  std::stringstream error_stream;
1656  error_stream
1657  << "Recast for element " << i << " faild" << std::endl;
1658  throw OomphLibError(
1659  error_stream.str(),
1660  OOMPH_CURRENT_FUNCTION,
1661  OOMPH_EXCEPTION_LOCATION);
1662  }
1663 #endif
1664 
1665  fe_pt->output_paraview(file_out,nplot);
1666  }
1667 
1668  file_out
1669  << "</DataArray>\n"
1670  << "</Points>\n";
1671 
1672 
1673  // Cells
1674  //-------
1675 
1676  file_out
1677  << "<Cells>\n"
1678  << "<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";
1679 
1680  // Make counter for keeping track of all the local elements,
1681  // because Paraview requires global coordinates
1682  unsigned counter=0;
1683 
1684  // Write connectivity with the local elements
1685  for(unsigned i=0;i<number_of_elements;i++)
1686  {
1687  // Cast to FiniteElement and (in paranoid mode) check
1688  // if cast has failed.
1689  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1690 
1691 #ifdef PARANOID
1692  if (fe_pt==0)
1693  {
1694  std::stringstream error_stream;
1695  error_stream
1696  << "Recast for element " << i << " faild" << std::endl;
1697  throw OomphLibError(
1698  error_stream.str(),
1699  OOMPH_CURRENT_FUNCTION,
1700  OOMPH_EXCEPTION_LOCATION);
1701  }
1702 #endif
1703  fe_pt->write_paraview_output_offset_information(file_out,nplot,counter);
1704  }
1705 
1706  file_out << "</DataArray>\n"
1707  << "<DataArray type=\"Int32\" "
1708  << "Name=\"offsets\" format=\"ascii\">\n";
1709 
1710  // Make variable that holds the current offset number
1711  unsigned offset_sum=0;
1712 
1713  // Write the offset for the specific elements
1714  for(unsigned i=0;i<number_of_elements;i++)
1715  {
1716  // Cast to FiniteElement and (in paranoid mode) check
1717  // if cast has failed.
1718  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1719 
1720 #ifdef PARANOID
1721  if (fe_pt==0)
1722  {
1723  std::stringstream error_stream;
1724  error_stream
1725  << "Recast for element " << i << " failed" << std::endl;
1726  throw OomphLibError(
1727  error_stream.str(),
1728  OOMPH_CURRENT_FUNCTION,
1729  OOMPH_EXCEPTION_LOCATION);
1730  }
1731 #endif
1732  fe_pt->write_paraview_offsets(file_out,nplot,offset_sum);
1733  }
1734 
1735  file_out <<"</DataArray>\n"
1736  <<"<DataArray type=\"UInt8\" Name=\"types\">\n";
1737 
1738  // Loop over all elements to get the type that they have
1739  for(unsigned i=0;i<number_of_elements;i++)
1740  {
1741  // Cast to FiniteElement and (in paranoid mode) check
1742  // if cast has failed.
1743  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(element_pt(i));
1744 
1745 #ifdef PARANOID
1746  if (fe_pt==0)
1747  {
1748  std::stringstream error_stream;
1749  error_stream
1750  << "Recast for element " << i << " failed" << std::endl;
1751  throw OomphLibError(
1752  error_stream.str(),
1753  OOMPH_CURRENT_FUNCTION,
1754  OOMPH_EXCEPTION_LOCATION);
1755  }
1756 #endif
1757 
1758  fe_pt->write_paraview_type(file_out,nplot);
1759  }
1760 
1761  file_out <<"</DataArray>\n"
1762  <<"</Cells>\n";
1763 
1764 
1765  // File Closure
1766  //-------------
1767  file_out <<"</Piece>\n"
1768  <<"</UnstructuredGrid>\n"
1769  <<"</VTKFile>";
1770 }
1771 
1772 //========================================================
1773 /// Output function for the mesh class
1774 ///
1775 /// Loop over all elements and plot (i.e. execute
1776 /// the element's own output() function)
1777 //========================================================
1778 void Mesh::output(std::ostream &outfile)
1779 {
1780  //Loop over the elements and call their output functions
1781  //Assign Element_pt_range
1782  unsigned long Element_pt_range = Element_pt.size();
1783  for(unsigned long e=0;e<Element_pt_range;e++)
1784  {
1785  // Try to cast to FiniteElement
1786  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
1787  if (el_pt==0)
1788  {
1789  oomph_info << "Can't execute output(...) for non FiniteElements"
1790  << std::endl;
1791  }
1792  else
1793  {
1794 #ifdef OOMPH_HAS_MPI
1796 #endif
1797  {
1798  el_pt->output(outfile);
1799  }
1800 #ifdef OOMPH_HAS_MPI
1801  else
1802  {
1803  if (!el_pt->is_halo())
1804  {
1805  el_pt->output(outfile);
1806  }
1807  }
1808 #endif
1809  }
1810  }
1811 }
1812 
1813 //========================================================
1814 /// Output function for the mesh class
1815 ///
1816 /// Loop over all elements and plot (i.e. execute
1817 /// the element's own output() function). Use
1818 /// n_plot plot points in each coordinate direction.
1819 //========================================================
1820 void Mesh::output(std::ostream &outfile, const unsigned &n_plot)
1821 {
1822  //Loop over the elements and call their output functions
1823  //Assign Element_pt_range
1824  unsigned long Element_pt_range = Element_pt.size();
1825 
1826  for(unsigned long e=0;e<Element_pt_range;e++)
1827  {
1828  // Try to cast to FiniteElement
1829  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
1830  if (el_pt==0)
1831  {
1832  oomph_info << "Can't execute output(...) for non FiniteElements"
1833  << std::endl;
1834  }
1835  else
1836  {
1837 #ifdef OOMPH_HAS_MPI
1839 #endif
1840  {
1841  el_pt->output(outfile,n_plot);
1842  }
1843 #ifdef OOMPH_HAS_MPI
1844  else
1845  {
1846  if (!el_pt->is_halo())
1847  {
1848  el_pt->output(outfile,n_plot);
1849  }
1850  }
1851 #endif
1852  }
1853  }
1854 }
1855 
1856 
1857 
1858 
1859 
1860 //========================================================
1861 /// Output function for the mesh class
1862 ///
1863 /// Loop over all elements and plot (i.e. execute
1864 /// the element's own output() function)
1865 /// (C style output)
1866 //========================================================
1867 void Mesh::output(FILE* file_pt)
1868 {
1869  //Loop over the elements and call their output functions
1870  //Assign Element_pt_range
1871  unsigned long Element_pt_range = Element_pt.size();
1872  for(unsigned long e=0;e<Element_pt_range;e++)
1873  {
1874  // Try to cast to FiniteElement
1875  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
1876  if (el_pt==0)
1877  {
1878  oomph_info << "Can't execute output(...) for non FiniteElements"
1879  << std::endl;
1880  }
1881  else
1882  {
1883 #ifdef OOMPH_HAS_MPI
1885 #endif
1886  {
1887  el_pt->output(file_pt);
1888  }
1889 #ifdef OOMPH_HAS_MPI
1890  else
1891  {
1892  if (!el_pt->is_halo())
1893  {
1894  el_pt->output(file_pt);
1895  }
1896  }
1897 #endif
1898  }
1899  }
1900 }
1901 
1902 //========================================================
1903 /// Output function for the mesh class
1904 ///
1905 /// Loop over all elements and plot (i.e. execute
1906 /// the element's own output() function). Use
1907 /// n_plot plot points in each coordinate direction.
1908 /// (C style output)
1909 //========================================================
1910 void Mesh::output(FILE* file_pt, const unsigned &n_plot)
1911 {
1912  //Loop over the elements and call their output functions
1913  //Assign Element_pt_range
1914  unsigned long Element_pt_range = Element_pt.size();
1915  for(unsigned long e=0;e<Element_pt_range;e++)
1916  {
1917  // Try to cast to FiniteElement
1918  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
1919  if (el_pt==0)
1920  {
1921  oomph_info << "Can't execute output(...) for non FiniteElements"
1922  << std::endl;
1923  }
1924  else
1925  {
1926 #ifdef OOMPH_HAS_MPI
1928 #endif
1929  {
1930  el_pt->output(file_pt,n_plot);
1931  }
1932 #ifdef OOMPH_HAS_MPI
1933  else
1934  {
1935  if (!el_pt->is_halo())
1936  {
1937  el_pt->output(file_pt,n_plot);
1938  }
1939  }
1940 #endif
1941  }
1942  }
1943 }
1944 
1945 
1946 //========================================================
1947 /// Output function for the mesh class
1948 ///
1949 /// Loop over all elements and plot (i.e. execute
1950 /// the element's own output() function). Use
1951 /// n_plot plot points in each coordinate direction.
1952 //========================================================
1953 void Mesh::output_fct(std::ostream &outfile, const unsigned &n_plot,
1955 {
1956  //Loop over the elements and call their output functions
1957  //Assign Element_pt_range
1958  unsigned long Element_pt_range = Element_pt.size();
1959  for(unsigned long e=0;e<Element_pt_range;e++)
1960  {
1961  // Try to cast to FiniteElement
1962  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
1963  if (el_pt==0)
1964  {
1965  oomph_info << "Can't execute output_fct(...) for non FiniteElements"
1966  << std::endl;
1967  }
1968  else
1969  {
1970 #ifdef OOMPH_HAS_MPI
1972 #endif
1973  {
1974  el_pt->output_fct(outfile,n_plot,exact_soln_pt);
1975  }
1976 #ifdef OOMPH_HAS_MPI
1977  else
1978  {
1979  if (!el_pt->is_halo())
1980  {
1981  el_pt->output_fct(outfile,n_plot,exact_soln_pt);
1982  }
1983  }
1984 #endif
1985  }
1986  }
1987 }
1988 
1989 //========================================================
1990 /// Output function for the mesh class
1991 ///
1992 /// Loop over all elements and plot (i.e. execute
1993 /// the element's own output() function) at time t. Use
1994 /// n_plot plot points in each coordinate direction.
1995 //========================================================
1996 void Mesh::output_fct(std::ostream &outfile, const unsigned &n_plot,
1997  const double& time,
1999 {
2000  //Loop over the elements and call their output functions
2001  //Assign Element_pt_range
2002  unsigned long Element_pt_range = Element_pt.size();
2003  for(unsigned long e=0;e<Element_pt_range;e++)
2004  {
2005  // Try to cast to FiniteElement
2006  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
2007  if (el_pt==0)
2008  {
2009  oomph_info << "Can't execute output_fct(...) for non FiniteElements"
2010  << std::endl;
2011  }
2012  else
2013  {
2014 #ifdef OOMPH_HAS_MPI
2016 #endif
2017  {
2018  el_pt->output_fct(outfile,n_plot,time,exact_soln_pt);
2019  }
2020 #ifdef OOMPH_HAS_MPI
2021  else
2022  {
2023  if (!el_pt->is_halo())
2024  {
2025  el_pt->output_fct(outfile,n_plot,time,exact_soln_pt);
2026  }
2027  }
2028 #endif
2029  }
2030  }
2031 }
2032 
2033 //==================================================================
2034 /// Assign the initial values for an impulsive start, which is
2035 /// acheived by looping over all data in the mesh (internal element
2036 /// data and data stored at nodes) and setting the calling the
2037 /// assign_initial_values_impulsive() function for each data's
2038 /// timestepper
2039 //=================================================================
2041 {
2042  //Loop over the elements
2043  unsigned long Nelement=nelement();
2044  for(unsigned long e=0;e<Nelement;e++)
2045  {
2046  //Find the number of internal dofs
2047  unsigned Ninternal = element_pt(e)->ninternal_data();
2048  //Loop over internal dofs and shift the time values
2049  //using the internals data's timestepper
2050  for(unsigned j=0;j<Ninternal;j++)
2051  {
2052  element_pt(e)->internal_data_pt(j)->time_stepper_pt()->
2053  assign_initial_values_impulsive(element_pt(e)->internal_data_pt(j));
2054  }
2055  }
2056 
2057  //Loop over the nodes
2058  unsigned long n_node=nnode();
2059  for (unsigned long n=0;n<n_node;n++)
2060  {
2061  // Assign initial values using the Node's timestepper
2062  Node_pt[n]->time_stepper_pt()->
2064  // Assign initial positions using the Node's timestepper
2065  Node_pt[n]->position_time_stepper_pt()->
2066  assign_initial_positions_impulsive(Node_pt[n]);
2067  }
2068 
2069 }
2070 
2071 //===============================================================
2072 /// Shift time-dependent data along for next timestep:
2073 /// Again this is achieved by looping over all data and calling
2074 /// the functions defined in each data object's timestepper.
2075 //==============================================================
2077 {
2078  // Loop over the elements which shift their internal data
2079  // via their own timesteppers
2080  const unsigned long Nelement=nelement();
2081  for (unsigned long e=0;e<Nelement;e++)
2082  {
2083  //Find the number of internal dofs
2084  const unsigned Ninternal = element_pt(e)->ninternal_data();
2085  //Loop over internal dofs and shift the time values
2086  //using the internals data's timestepper
2087  for(unsigned j=0;j<Ninternal;j++)
2088  {
2089  element_pt(e)->internal_data_pt(j)->time_stepper_pt()->
2090  shift_time_values(element_pt(e)->internal_data_pt(j));
2091  }
2092  }
2093 
2094  //Loop over the nodes
2095  const unsigned long n_node=nnode();
2096  for (unsigned long n=0;n<n_node;n++)
2097  {
2098  // Shift the Data associated with the nodes with the Node's own
2099  // timestepper
2100  Node_pt[n]->time_stepper_pt()->shift_time_values(Node_pt[n]);
2101  // Push history of nodal positions back
2102  Node_pt[n]->position_time_stepper_pt()->shift_time_positions(Node_pt[n]);
2103  }
2104 
2105 }
2106 
2107 //=========================================================================
2108 /// Calculate predictions for all Data and positions associated
2109 /// with the mesh. This is usually only used for adaptive time-stepping
2110 /// when the comparison between a predicted value and the actual value
2111 /// is usually used to determine the change in step size. Again the
2112 /// loop is over all data in the mesh and individual timestepper functions
2113 /// for each data value are called.
2114 //=========================================================================
2116 {
2117  // Loop over the elements which shift their internal data
2118  // via their own timesteppers
2119  const unsigned long Nelement=nelement();
2120  for (unsigned long e=0;e<Nelement;e++)
2121  {
2122  //Find the number of internal dofs
2123  const unsigned Ninternal = element_pt(e)->ninternal_data();
2124  //Loop over internal dofs and calculate predicted positions
2125  //using the internals data's timestepper
2126  for(unsigned j=0;j<Ninternal;j++)
2127  {
2128  element_pt(e)->internal_data_pt(j)->time_stepper_pt()->
2129  calculate_predicted_values(element_pt(e)->internal_data_pt(j));
2130  }
2131  }
2132 
2133  //Loop over the nodes
2134  const unsigned long n_node=nnode();
2135  for (unsigned long n=0;n<n_node;n++)
2136  {
2137  // Calculate the predicted values at the nodes
2138  Node_pt[n]->time_stepper_pt()->calculate_predicted_values(Node_pt[n]);
2139  //Calculate the predicted positions
2140  Node_pt[n]->position_time_stepper_pt()->
2141  calculate_predicted_positions(Node_pt[n]);
2142  }
2143 }
2144 
2145 //===============================================================
2146 /// Virtual function that should be overloaded if the mesh
2147 /// has any mesh level storage of the timestepper
2148 //==================================================================
2149 void Mesh::set_mesh_level_time_stepper(TimeStepper* const &time_stepper_pt,
2150  const bool &preserve_existing_data)
2151 {
2152 #ifdef PARANOID
2154  {
2155  std::ostringstream warning_stream;
2156  warning_stream <<
2157  "Empty set_mesh_level_time_stepper() has been called.\n"
2158  <<
2159  "This function needs to be overloaded to reset any (pointers to) \n"
2160  <<
2161  "timesteppers for meshes that store timesteppers in locations other\n"
2162  << "than the Nodes or Elements;\n"
2163  << "e.g. SpineMeshes have SpineData with timesteppers,\n"
2164  <<
2165  "Triangle and TetMeshes store the timestepper for use in adaptivity.\n\n\n";
2166  warning_stream << "If you are solving a continuation or bifurcation detecion\n"
2167  << "problem and strange things are happening, then check that\n"
2168  << "you don't need to overload this function for your mesh."
2169  << "\n This warning can be suppressed by setting:\n"
2170  <<
2171  "Mesh::Suppress_warning_about_empty_mesh_level_time_stepper_function=true"
2172  << std::endl;
2173  OomphLibWarning(warning_stream.str(),
2174  OOMPH_CURRENT_FUNCTION,
2175  OOMPH_EXCEPTION_LOCATION);
2176  }
2177 #endif
2178 }
2179 
2180 //===============================================================
2181 /// Set the values of auxilliary data used in continuation problems
2182 /// when the data is pinned.
2183 //==================================================================
2185  ContinuationStorageScheme* const &continuation_storage_pt)
2186 {
2187  //Loop over the nodes
2188  const unsigned long n_node=this->nnode();
2189  for(unsigned long n=0;n<n_node;n++)
2190  {
2191  continuation_storage_pt->set_consistent_pinned_values(this->Node_pt[n]);
2192  continuation_storage_pt->set_consistent_pinned_positions(this->Node_pt[n]);
2193  }
2194 
2195  // Loop over the elements
2196  const unsigned long n_element=this->nelement();
2197  for (unsigned long e=0;e<n_element;e++)
2198  {
2199  //Cache pointer to the elemnet
2200  GeneralisedElement* const elem_pt = this->element_pt(e);
2201  //Find the number of internal dofs
2202  const unsigned n_internal = elem_pt->ninternal_data();
2203 
2204  //Loop over internal dofs and test the data
2205  for(unsigned j=0;j<n_internal;j++)
2206  {
2207  continuation_storage_pt->
2208  set_consistent_pinned_values(elem_pt->internal_data_pt(j));
2209  }
2210  }
2211 }
2212 
2213 
2214 //===============================================================
2215 /// Return true if the pointer corresponds to any data stored in
2216 /// the mesh and false if not
2217 //==================================================================
2218 bool Mesh::does_pointer_correspond_to_mesh_data(double* const &parameter_pt)
2219 {
2220  //Loop over the nodes
2221  const unsigned long n_node=this->nnode();
2222  for(unsigned long n=0;n<n_node;n++)
2223  {
2224  //Check the values and positional data associated with each node
2225  if(
2226  (this->Node_pt[n]->does_pointer_correspond_to_value(parameter_pt))
2227  ||
2228  (this->Node_pt[n]->does_pointer_correspond_to_position_data(parameter_pt)))
2229  {return true;}
2230  }
2231 
2232  // Loop over the elements
2233  const unsigned long n_element=this->nelement();
2234  for (unsigned long e=0;e<n_element;e++)
2235  {
2236  //Cache pointer to the elemnet
2237  GeneralisedElement* const elem_pt = this->element_pt(e);
2238 
2239  //Find the number of internal dofs
2240  const unsigned n_internal = elem_pt->ninternal_data();
2241 
2242  //Loop over internal dofs and test the data
2243  for(unsigned j=0;j<n_internal;j++)
2244  {
2245  if(elem_pt->internal_data_pt(j)
2246  ->does_pointer_correspond_to_value(parameter_pt))
2247  {return true;}
2248  }
2249  }
2250 
2251  //If we get here we haven't found the data, so return false
2252  return false;
2253 }
2254 
2255 
2256 
2257 //===============================================================
2258 /// Set the time stepper associated with all the nodal data
2259 /// in the problem
2260 //==============================================================
2261 void Mesh::set_nodal_time_stepper(TimeStepper* const &time_stepper_pt,
2262  const bool &preserve_existing_data)
2263 {
2264  //Loop over the nodes
2265  const unsigned long n_node=this->nnode();
2266  for(unsigned long n=0;n<n_node;n++)
2267  {
2268  //Set the timestepper associated with each node
2269  this->Node_pt[n]->set_time_stepper(time_stepper_pt,preserve_existing_data);
2270  this->Node_pt[n]->set_position_time_stepper(time_stepper_pt,
2271  preserve_existing_data);
2272  }
2273 }
2274 
2275 //===============================================================
2276 /// Set the time stepper associated with all internal data stored
2277 /// in the elements in the mesh
2278 //===============================================================
2280  TimeStepper* const &time_stepper_pt, const bool &preserve_existing_data)
2281 {
2282  // Loop over the elements
2283  const unsigned long n_element=this->nelement();
2284  for (unsigned long e=0;e<n_element;e++)
2285  {
2286  //Find the number of internal dofs
2287  const unsigned n_internal = this->element_pt(e)->ninternal_data();
2288 
2289  //Loop over internal dofs and set the timestepper
2290  for(unsigned j=0;j<n_internal;j++)
2291  {
2292  this->element_pt(e)->internal_data_pt(j)
2293  ->set_time_stepper(time_stepper_pt,preserve_existing_data);
2294  }
2295  }
2296 }
2297 
2298 //========================================================================
2299 /// A function that upgrades an ordinary node to a boundary node.
2300 /// All pointers to the node from the mesh's elements are found.
2301 /// and replaced by pointers to the new boundary node. If the node
2302 /// is present in the mesh's list of nodes, that pointer is also
2303 /// replaced. Finally, the pointer argument node_pt addresses the new
2304 /// node on return from the function.
2305 /// We shouldn't ever really use this, but it does make life that
2306 /// bit easier for the lazy mesh writer.
2307 //=======================================================================
2309 {
2310  // Cache a list of FiniteElement pointers for use in this function.
2311  Vector<FiniteElement*> fe_pt(nelement(),0);
2312  for(unsigned e=0, ne=nelement(); e < ne; e++)
2313  {
2314  // Some elements may not have been build yet, just store a null pointer
2315  // for these cases.
2316  if(Element_pt[e] == 0) fe_pt[e] = 0;
2317  else fe_pt[e] = finite_element_pt(e);
2318  }
2319 
2320  // Now call the real function
2321  convert_to_boundary_node(node_pt, fe_pt);
2322 }
2323 
2324 // ============================================================
2325 /// As convert_to_boundary_node but with a vector of pre-"dynamic cast"ed
2326 /// pointers passed in. If this function is being called often then
2327 /// creating this vector and passing it in explicitly each time can give a
2328 /// large speed up.
2329 // Note: the real reason that this function is so slow in the first place
2330 // is because it has to loop over all elements. So if you use this function
2331 // O(N) times your boundary node creation complexity is O(N^2).
2332 // ============================================================
2334 (Node* &node_pt, const Vector<FiniteElement*>& finite_element_pt)
2335 {
2336 
2337  //If the node is already a boundary node, then return straight away,
2338  //we don't need to do anything
2339  if (dynamic_cast<BoundaryNodeBase*>(node_pt)!=0)
2340  {
2341  return;
2342  }
2343 
2344  //Loop over all the elements in the mesh and find all those in which
2345  //the present node is referenced and the corresponding local node number
2346  //in those elements.
2347 
2348  //Storage for elements and local node number
2349  std::list<std::pair<unsigned long, int> >
2350  list_of_elements_and_local_node_numbers;
2351 
2352  //Loop over all elements
2353  unsigned long n_element=this->nelement();
2354  for(unsigned long e=0;e<n_element;e++)
2355  {
2356  //Buffer the case when we have not yet filled up the element array
2357  //Unfortunately, we should not assume that the array has been filled
2358  //in a linear order, so we can't break out early.
2359  if(Element_pt[e]!=0)
2360  {
2361  //Find the local node number of the passed node
2362  int node_number = finite_element_pt[e]->get_node_number(node_pt);
2363  //If the node is present in the element, add it to our list and
2364  //NULL out the local element entries
2365  if(node_number!=-1)
2366  {
2367  list_of_elements_and_local_node_numbers.insert(
2368  list_of_elements_and_local_node_numbers.end(),
2369  std::make_pair(e,node_number));
2370  //Null it out
2371  finite_element_pt[e]->node_pt(node_number)=0;
2372  }
2373  }
2374  } //End of loop over elements
2375 
2376  //If there are no entries in the list we are in real trouble
2377  if(list_of_elements_and_local_node_numbers.empty())
2378  {
2379  std::ostringstream error_stream;
2380  error_stream << "Node " << node_pt
2381  << " is not contained in any elements in the Mesh."
2382  << std::endl
2383  << "How was it created then?" << std::endl;
2384 
2385  throw OomphLibError(error_stream.str(),
2386  OOMPH_CURRENT_FUNCTION,
2387  OOMPH_EXCEPTION_LOCATION);
2388  }
2389 
2390 
2391  //Create temporary storage for a pointer to the old node.
2392  //This is required because if we have passed a reference to the
2393  //first element that we find, constructing the new node
2394  //will over-write our pointer and we'll get segmentation faults.
2395  Node* old_node_pt = node_pt;
2396 
2397  //We now create the new node by using the first element in the list
2398  std::list<std::pair<unsigned long,int> >::iterator list_it
2399  = list_of_elements_and_local_node_numbers.begin();
2400 
2401  //Create a new boundary node, using the timestepper from the
2402  //original node
2403  Node* new_node_pt = finite_element_pt[list_it->first]
2404  ->construct_boundary_node(list_it->second,node_pt->time_stepper_pt());
2405 
2406  //Now copy all the information accross from the old node
2407 
2408  //Can we cast the node to a solid node
2409  SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(new_node_pt);
2410  //If it's a solid node, do the casting
2411  if(solid_node_pt!=0)
2412  {
2413  solid_node_pt->copy(dynamic_cast<SolidNode*>(old_node_pt));
2414  }
2415  else
2416  {
2417  new_node_pt->copy(old_node_pt);
2418  }
2419 
2420  //Loop over all other elements in the list and set their pointers
2421  //to the new node
2422  for(++list_it; //Increment the iterator
2423  list_it!=list_of_elements_and_local_node_numbers.end();
2424  ++list_it)
2425  {
2426  finite_element_pt[list_it->first]->node_pt(list_it->second)
2427  = new_node_pt;
2428  }
2429 
2430  //Finally, find the position of the node in the global mesh
2432  std::find(Node_pt.begin(),Node_pt.end(),old_node_pt);
2433 
2434  //If it is in the mesh, update the pointer
2435  if(it!=Node_pt.end()) {*it = new_node_pt;}
2436 
2437  //Can now delete the old node
2438  delete old_node_pt;
2439 
2440  //Replace the passed pointer by a pointer to the new node
2441  //Note that in most cases, this will be wasted work because the node
2442  //pointer will either the pointer in the mesh's or an element's
2443  //node_pt vector. Still assignment is quicker than an if to check this.
2444  node_pt = new_node_pt;
2445 
2446 }
2447 
2448 #ifdef OOMPH_HAS_MPI
2449 
2450 //========================================================================
2451 /// Setup shared node scheme: Shared node lookup scheme contains
2452 /// a unique correspondence between all nodes on the halo/haloed
2453 /// elements between two processors.
2454 //========================================================================
2456 {
2457 
2458  // Determine the shared nodes lookup scheme - all nodes located on the
2459  // halo(ed) elements between two domains. This scheme is necessary in order
2460  // to identify master nodes that may not be present in the halo-haloed
2461  // element lookup scheme between two processors (for example, if the node
2462  // is on an element which is in a lookup scheme between two higher-numbered
2463  // processors)
2464 
2465  double t_start = 0.0;
2467  {
2468  t_start=TimingHelpers::timer();
2469  }
2470 
2471  // Store number of processors and current process
2472  int n_proc=Comm_pt->nproc();
2473  int my_rank=Comm_pt->my_rank();
2474 
2475  // Need to clear the shared node scheme first
2476  Shared_node_pt.clear();
2477 
2478  for (int d=0;d<n_proc;d++)
2479  {
2480 
2481  // map of bools for whether the node has been shared,
2482  // initialised to 0 (false) for each domain d
2483  std::map<Node*,bool> node_shared;
2484 
2485  // For all domains lower than the current domain: Do halos first
2486  // then haloed, to ensure correct order in lookup scheme from
2487  // the other side
2488  if (d<my_rank)
2489  {
2490  // Get the nodes from the halo elements first
2491  Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(d));
2492  unsigned nhalo_elem=halo_elem_pt.size();
2493 
2494  for (unsigned e=0;e<nhalo_elem;e++)
2495  {
2496  // Get finite element
2497  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(halo_elem_pt[e]);
2498  if (el_pt!=0)
2499  {
2500  // Loop over nodes
2501  unsigned nnod=el_pt->nnode();
2502  for (unsigned j=0;j<nnod;j++)
2503  {
2504  Node* nod_pt=el_pt->node_pt(j);
2505 
2506  // Add it as a shared node from current domain
2507  if (!node_shared[nod_pt])
2508  {
2509  this->add_shared_node_pt(d,nod_pt);
2510  node_shared[nod_pt]=true;
2511  }
2512 
2513  } // end loop over nodes
2514  }
2515 
2516  } // end loop over elements
2517 
2518  // Now get any left over nodes on the haloed elements
2519  Vector<GeneralisedElement*> haloed_elem_pt(this->haloed_element_pt(d));
2520  unsigned nhaloed_elem=haloed_elem_pt.size();
2521 
2522  for (unsigned e=0;e<nhaloed_elem;e++)
2523  {
2524  // Get element
2525  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(haloed_elem_pt[e]);
2526  if (el_pt!=0)
2527  {
2528  // Loop over the nodes
2529  unsigned nnod=el_pt->nnode();
2530  for (unsigned j=0;j<nnod;j++)
2531  {
2532  Node* nod_pt=el_pt->node_pt(j);
2533 
2534  // Add it as a shared node from current domain
2535  if (!node_shared[nod_pt])
2536  {
2537  this->add_shared_node_pt(d,nod_pt);
2538  node_shared[nod_pt]=true;
2539  }
2540 
2541  } // end loop over nodes
2542  }
2543  } // end loop over elements
2544 
2545  }
2546 
2547  // If the domain is bigger than the current rank: Do haloed first
2548  // then halo, to ensure correct order in lookup scheme from
2549  // the other side
2550  if (d>my_rank)
2551  {
2552  // Get the nodes from the haloed elements first
2553  Vector<GeneralisedElement*> haloed_elem_pt(this->haloed_element_pt(d));
2554  unsigned nhaloed_elem=haloed_elem_pt.size();
2555 
2556  for (unsigned e=0;e<nhaloed_elem;e++)
2557  {
2558  // Get element
2559  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(haloed_elem_pt[e]);
2560  if (el_pt!=0)
2561  {
2562  // Loop over nodes
2563  unsigned nnod=el_pt->nnode();
2564  for (unsigned j=0;j<nnod;j++)
2565  {
2566  Node* nod_pt=el_pt->node_pt(j);
2567 
2568  // Add it as a shared node from current domain
2569  if (!node_shared[nod_pt])
2570  {
2571  this->add_shared_node_pt(d,nod_pt);
2572  node_shared[nod_pt]=true;
2573  }
2574 
2575  } // end loop over nodes
2576  }
2577  } // end loop over elements
2578 
2579  // Now get the nodes from any halo elements left over
2580  Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(d));
2581  unsigned nhalo_elem=halo_elem_pt.size();
2582 
2583  for (unsigned e=0;e<nhalo_elem;e++)
2584  {
2585  // Get element
2586  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(halo_elem_pt[e]);
2587  if (el_pt!=0)
2588  {
2589  // Loop over nodes
2590  unsigned nnod=el_pt->nnode();
2591  for (unsigned j=0;j<nnod;j++)
2592  {
2593  Node* nod_pt=el_pt->node_pt(j);
2594 
2595  // Add it as a shared node from current domain
2596  if (!node_shared[nod_pt])
2597  {
2598  this->add_shared_node_pt(d,nod_pt);
2599  node_shared[nod_pt]=true;
2600  }
2601 
2602  } // end loop over nodes
2603  }
2604  } // end loop over elements
2605 
2606  } // end if (d ...)
2607 
2608  } // end loop over processes
2609 
2610 
2611  double t_end=0.0;
2613  {
2614  t_end = TimingHelpers::timer();
2615  oomph_info << "Time for identification of shared nodes: "
2616  << t_end-t_start << std::endl;
2617  oomph_info.stream_pt()->flush();
2618  }
2619 
2620 }
2621 
2622 //========================================================================
2623 /// \short Synchronise shared node lookup schemes to cater for the
2624 /// the case where:
2625 /// (1) a certain node on the current processor is halo with proc p
2626 /// (i.e. its non-halo counterpart lives on processor p)
2627 /// (2) that node also exists (also as a halo) on another processor
2628 /// (q, say) where its non-halo counter part is also known to be
2629 /// on processor p.
2630 /// However, without calling this function the current processor does not
2631 /// necessarily know that it shares a node with processor q. This
2632 /// information can be required, e.g. when synchronising hanging node
2633 /// schemes over all processors.
2634 //========================================================================
2635 void Mesh::synchronise_shared_nodes(const bool& report_stats)
2636 {
2637 
2638  double t_start=0.0;
2640  {
2641  t_start=TimingHelpers::timer();
2642  }
2643 
2644  double tt_start=0.0;
2645  double tt_end=0.0;
2647  {
2648  tt_start=TimingHelpers::timer();
2649  }
2650 
2651  // Storage for current processor and number of processors
2652  int n_proc=Comm_pt->nproc();
2653  int my_rank=Comm_pt->my_rank();
2654 
2655 
2656 #ifdef PARANOID
2657  // Has some twit filled in shared nodes with own process?!
2658  // Check at start of function
2659  if (Shared_node_pt[my_rank].size()!=0)
2660  {
2661  throw OomphLibError(
2662  "Processor has shared nodes with itself! Something's gone wrong!",
2663  OOMPH_CURRENT_FUNCTION,
2664  OOMPH_EXCEPTION_LOCATION);
2665  }
2666 #endif
2667 
2668 
2669  // Stage 1: Populate the set of of processor IDs that
2670  // each haloed node on current processor is haloed by.
2671  std::map<Node*,std::set<int> > shared_domain_set;
2672 
2673  // Associate unique number with any haloed nodes on this processor
2674  std::map<Node*,unsigned> global_haloed_node_number_plus_one;
2675  unsigned global_haloed_count=0;
2676 
2677  // Loop over domains
2678  for (int d=0;d<n_proc;d++)
2679  {
2680  // Don't talk to yourself
2681  if (d!=my_rank)
2682  {
2683  // Loop over haloed nodes
2684  unsigned nnod_haloed=this->nhaloed_node(d);
2685  for (unsigned j=0;j<nnod_haloed;j++)
2686  {
2687  Node* nod_pt=this->haloed_node_pt(d,j);
2688  shared_domain_set[nod_pt].insert(d);
2689  if (global_haloed_node_number_plus_one[nod_pt]==0)
2690  {
2691  global_haloed_node_number_plus_one[nod_pt]=global_haloed_count+1;
2692  global_haloed_count++;
2693  }
2694  }
2695  }
2696  }
2697 
2699  {
2700  tt_end = TimingHelpers::timer();
2701  oomph_info
2702  << "Time for initial classification in Mesh::synchronise_shared_nodes(): "
2703  << tt_end-tt_start << std::endl;
2704  tt_start = TimingHelpers::timer();
2705  }
2706 
2707 
2708  // Stage 2: All-to-all communication to inform all processors that
2709  // hold halo nodes with current processor about all domains that the current
2710  // processor shares nodes with [This will allow the other processors to add
2711  // these nodes to their shared node lookup schemes with processors
2712  // that they currently "can't see"].
2713 
2714  // Data to be sent to each processor
2715  Vector<int> send_n(n_proc,0);
2716 
2717  // Storage for all values to be sent to all processors
2718  Vector<unsigned> send_data;
2719 
2720  // Start location within send_data for data to be sent to each processor
2721  Vector<int> send_displacement(n_proc,0);
2722 
2723 
2724  // Loop over haloed nodes with other domains
2725  for (int domain=0;domain<n_proc;domain++)
2726  {
2727  //Set the offset for the current processor
2728  send_displacement[domain] = send_data.size();
2729 
2730  // Every processor works on haloed nodes with proc "domain" and
2731  // sends associations across to that domain. No need to talk to yourself...
2732  if(domain!=my_rank)
2733  {
2734  // Send total number of global haloed nodes
2735  // send_data.push_back(global_haloed_count);
2736 
2737  // Loop over haloed nodes
2738  unsigned nnod_haloed=this->nhaloed_node(domain);
2739  for (unsigned j=0;j<nnod_haloed;j++)
2740  {
2741  Node* nod_pt=this->haloed_node_pt(domain,j);
2742 
2743  // Send global ID of haloed node
2744  send_data.push_back(global_haloed_node_number_plus_one[nod_pt]-1);
2745 
2746  // Get the set of domains that halo this node
2747  std::set<int> tmp_shared_domain_set=shared_domain_set[nod_pt];
2748 
2749  // Insert number of shared domains into send data
2750  unsigned n_shared_domains= tmp_shared_domain_set.size();
2751  send_data.push_back(n_shared_domains);
2752 
2753  // Add shared domains
2754  for (std::set<int>::iterator it=tmp_shared_domain_set.begin();
2755  it!=tmp_shared_domain_set.end();it++)
2756  {
2757  send_data.push_back(*it);
2758  }
2759  }
2760  }
2761 
2762  //Find the number of data added to the vector
2763  send_n[domain] = send_data.size() - send_displacement[domain];
2764 
2765  }
2766 
2767 
2768 
2769  //Storage for the number of data to be received from each processor
2770  Vector<int> receive_n(n_proc,0);
2771 
2772  //Now send numbers of data to be sent between all processors
2773  MPI_Alltoall(&send_n[0],1,MPI_INT,&receive_n[0],1,MPI_INT,
2774  Comm_pt->mpi_comm());
2775 
2776 
2777  //We now prepare the data to be received
2778  //by working out the displacements from the received data
2779  Vector<int> receive_displacement(n_proc,0);
2780  int receive_data_count=0;
2781  for(int rank=0;rank<n_proc;++rank)
2782  {
2783  //Displacement is number of data received so far
2784  receive_displacement[rank] = receive_data_count;
2785  receive_data_count += receive_n[rank];
2786  }
2787 
2788  //Now resize the receive buffer for all data from all processors
2789  //Make sure that it has a size of at least one
2790  if(receive_data_count==0) {++receive_data_count;}
2791  Vector<unsigned> receive_data(receive_data_count);
2792 
2793  //Make sure that the send buffer has size at least one
2794  //so that we don't get a segmentation fault
2795  if(send_data.size()==0) {send_data.resize(1);}
2796 
2797  //Now send the data between all the processors
2798  MPI_Alltoallv(&send_data[0],&send_n[0],&send_displacement[0],
2799  MPI_UNSIGNED,
2800  &receive_data[0],&receive_n[0],
2801  &receive_displacement[0],
2802  MPI_UNSIGNED,
2803  Comm_pt->mpi_comm());
2804 
2805 
2807  {
2808  tt_end = TimingHelpers::timer();
2809  oomph_info
2810  << "Time for alltoall in Mesh::synchronise_shared_nodes(): "
2811  << tt_end-tt_start << std::endl;
2812  tt_start = TimingHelpers::timer();
2813  }
2814 
2815 
2817  {
2818  oomph_info
2819  << "Starting vector to set conversion in "
2820  << "Mesh::synchronise_shared_nodes() for a total of "
2821  << nshared_node() << " nodes\n";
2822  tt_start = TimingHelpers::timer();
2823  }
2824 
2825 
2826  // Copy vector-based representation of shared nodes into
2827  // sets for faster search
2828  Vector<std::set<Node*> > shared_node_set(n_proc);
2829  for (int d=0;d<n_proc;d++)
2830  {
2831  unsigned n_vector=Shared_node_pt[d].size();
2832  for (unsigned i=0;i<n_vector;i++)
2833  {
2834  shared_node_set[d].insert(Shared_node_pt[d][i]);
2835  }
2836  }
2837 
2838 
2840  {
2841  tt_end = TimingHelpers::timer();
2842  oomph_info
2843  << "Time for vector to set in Mesh::synchronise_shared_nodes(): "
2844  << tt_end-tt_start << std::endl;
2845  tt_start = TimingHelpers::timer();
2846  }
2847 
2848 
2849  //Now use the received data
2850  for (int send_rank=0;send_rank<n_proc;send_rank++)
2851  {
2852  //Don't bother to do anything for the processor corresponding to the
2853  //current processor or if no data were received from this processor
2854  if((send_rank != my_rank) && (receive_n[send_rank] != 0))
2855  {
2856  //Counter for the data within the large array
2857  unsigned count=receive_displacement[send_rank];
2858 
2859  // Read total number of global haloed nodes
2860  //unsigned n_global_haloed_nodes_on_send_proc=receive_data[count++];
2861 
2862  // Storage for nodes and associated domains:
2863  // domains_map[global_haloed_node_number].first = node
2864  // domains_map[global_haloed_node_number].second = set of domains
2865  // this node is associated
2866  // with.
2867  std::map<unsigned,std::pair<Node*,std::set<unsigned> > > domains_map;
2868 
2869  // Loop over halo nodes with sending processor
2870  unsigned nnod_halo=this->nhalo_node(send_rank);
2871  for (unsigned j=0;j<nnod_halo;j++)
2872  {
2873  Node* nod_pt=this->halo_node_pt(send_rank,j);
2874 
2875  // Read unique ID of haloed node on send proc
2876  unsigned haloed_node_id_on_send_proc=receive_data[count++];
2877 
2878  // Read out number of shared domains into send data
2879  unsigned n_shared_domains=receive_data[count++];
2880 
2881  // Prepare set of domains
2882  std::set<unsigned> domain_set;
2883 
2884  // Read 'em
2885  for (unsigned i=0;i<n_shared_domains;i++)
2886  {
2887  int shared_domain=receive_data[count++];
2888 
2889  // Record in set
2890  domain_set.insert(shared_domain);
2891  }
2892 
2893  // Add entry:
2894  domains_map[haloed_node_id_on_send_proc]=std::make_pair(nod_pt,
2895  domain_set);
2896 
2897  } // end of loop over halo nodes
2898 
2899 
2900 
2901  // Now add new shared nodes in order
2902 #ifdef PARANOID
2903  int previous_one=-1;
2904 #endif
2905  for (std::map<unsigned,std::pair<Node*,std::set<unsigned> > >::iterator
2906  it=domains_map.begin();it!=domains_map.end();it++)
2907  {
2908 
2909  // Super-paranoid: Check that the map actually sorted entries
2910  // by key (as it should)
2911 #ifdef PARANOID
2912  if (int((*it).first)<previous_one)
2913  {
2914  std::ostringstream error_stream;
2915  error_stream << "Map did not store entries in order of key\n "
2916  << "Current key: " << (*it).first
2917  << "; previous one: " << previous_one << "\n"
2918  << "Need to rewrite this...\n";
2919  throw OomphLibError(error_stream.str(),
2920  OOMPH_CURRENT_FUNCTION,
2921  OOMPH_EXCEPTION_LOCATION);
2922  }
2923  previous_one=(*it).first;
2924 #endif
2925 
2926  // Extract node
2927  Node* nod_pt=(*it).second.first;
2928 
2929  // Extract set of domains
2930  std::set<unsigned> domain_set((*it).second.second);
2931 
2932  // Read 'em
2933  for (std::set<unsigned>::iterator itt=domain_set.begin();
2934  itt!=domain_set.end();itt++)
2935  {
2936  int shared_domain=(*itt);
2937 
2938  // No need to add shared nodes with oneself!
2939  if (shared_domain!=my_rank)
2940  {
2941  // Is node already listed in shared node scheme? Find it
2942  // and get iterator to entry
2943  std::set<Node*>::iterator ittt=
2944  shared_node_set[shared_domain].find(nod_pt);
2945 
2946  // If it's not in there already iterator points to end of
2947  // set
2948  if (ittt==shared_node_set[shared_domain].end())
2949  {
2950  // Now add it
2951  add_shared_node_pt(shared_domain,nod_pt);
2952 
2953  // Update set
2954  shared_node_set[shared_domain].insert(nod_pt);
2955  }
2956  }
2957  }
2958  }
2959 
2960  } // end of any data is received and ignore own domain
2961 
2962  } // end of loop over send ranks
2963 
2964 
2965 
2966 #ifdef PARANOID
2967  // Has some twit filled in shared nodes with own process?!
2968  // Check at end pf function.
2969  if (Shared_node_pt[my_rank].size()!=0)
2970  {
2971  throw OomphLibError(
2972  "Processor has shared nodes with itself! Something's gone wrong!",
2973  OOMPH_CURRENT_FUNCTION,
2974  OOMPH_EXCEPTION_LOCATION);
2975  }
2976 #endif
2977 
2978 
2980  {
2981  tt_end = TimingHelpers::timer();
2982  oomph_info
2983  << "Time for final processing in Mesh::synchronise_shared_nodes(): "
2984  << tt_end-tt_start << std::endl;
2985  tt_start = TimingHelpers::timer();
2986  }
2987 
2989  {
2990  double t_end = TimingHelpers::timer();
2991  oomph_info << "Total time for Mesh::synchronise_shared_nodes(): "
2992  << t_end-t_start << std::endl;
2993  }
2994 }
2995 
2996 
2997 
2998 
2999 
3000 
3001 //========================================================================
3002 /// Classify all halo and haloed information in the mesh
3003 //========================================================================
3005  const bool& report_stats)
3006 {
3007 
3008 
3009 // MemoryUsage::doc_memory_usage(
3010 // "at beginning of Mesh::classify_halo_and_haloed_nodes");
3011 
3012  double t_start = 0.0;
3014  {
3015  t_start=TimingHelpers::timer();
3016  }
3017 
3018  double tt_start = 0.0;
3020  {
3021  tt_start=TimingHelpers::timer();
3022  }
3023 
3024  // Set up shared nodes scheme
3026 
3027  //MemoryUsage::doc_memory_usage("after setup shared node scheme");
3028 
3029  double tt_end=0.0;
3031  {
3032  tt_end=TimingHelpers::timer();
3033  oomph_info
3034  << "Time for Mesh::setup_shared_node_scheme() "
3035  << " Mesh::classify_halo_and_haloed_nodes(): "
3036  << tt_end-tt_start << std::endl;
3037  oomph_info.stream_pt()->flush();
3038  tt_start = TimingHelpers::timer();
3039  }
3040 
3041 
3042  //Wipe existing storage schemes for halo(ed) nodes
3043  Halo_node_pt.clear();
3044  Haloed_node_pt.clear();
3045 
3046  // Storage for number of processors and current processor
3047  int n_proc=Comm_pt->nproc();
3048  int my_rank=Comm_pt->my_rank();
3049  MPI_Status status;
3050 
3051  // Determine which processors the nodes are associated with
3052  // and hence who's in charge
3053  std::map<Data*,std::set<unsigned> > processors_associated_with_data;
3054  std::map<Data*,unsigned> processor_in_charge;
3055 
3056  // Loop over all processors and associate any nodes on the halo
3057  // elements involved with that processor
3058  for (int domain=0;domain<n_proc;domain++)
3059  {
3060  // Get vector of halo elements by copy operation
3061  Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(domain));
3062 
3063  // Loop over halo elements associated with this adjacent domain
3064  unsigned nelem=halo_elem_pt.size();
3065  for (unsigned e=0;e<nelem;e++)
3066  {
3067  // Get element only have nodes if a finite element
3068  FiniteElement* finite_el_pt=dynamic_cast<FiniteElement*>(halo_elem_pt[e]);
3069  if(finite_el_pt!=0)
3070  {
3071  //Loop over nodes
3072  unsigned nnod=finite_el_pt->nnode();
3073  for (unsigned j=0;j<nnod;j++)
3074  {
3075  Node* nod_pt=finite_el_pt->node_pt(j);
3076  // Associate node with this domain
3077  processors_associated_with_data[nod_pt].insert(domain);
3078 
3079  // Do the same if the node is solid
3080  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
3081  if (solid_nod_pt!=0)
3082  {
3083  processors_associated_with_data[
3084  solid_nod_pt->variable_position_pt()].insert(domain);
3085  }
3086  }
3087  }
3088  }
3089  }
3090 
3091 
3092  // Loop over all [non-halo] elements and associate their nodes
3093  // with current procesor
3094  unsigned nelem=this->nelement();
3095  for (unsigned e=0;e<nelem;e++)
3096  {
3097  FiniteElement* finite_el_pt=
3098  dynamic_cast<FiniteElement*>(this->element_pt(e));
3099 
3100  // Only visit non-halos and finite elements
3101  if((finite_el_pt!=0) && (!finite_el_pt->is_halo()))
3102  {
3103  // Loop over nodes
3104  unsigned nnod=finite_el_pt->nnode();
3105  for (unsigned j=0;j<nnod;j++)
3106  {
3107  Node* nod_pt=finite_el_pt->node_pt(j);
3108 
3109  // Associate this node with current processor
3110  processors_associated_with_data[nod_pt].insert(my_rank);
3111 
3112  // do the same if we have a SolidNode
3113  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
3114  if (solid_nod_pt!=0)
3115  {
3116  processors_associated_with_data
3117  [solid_nod_pt->variable_position_pt()].insert(my_rank);
3118  }
3119  }
3120  }
3121  }
3122 
3123 
3125  {
3126  tt_end = TimingHelpers::timer();
3127  oomph_info
3128  << "Time for setup loops in Mesh::classify_halo_and_haloed_nodes: "
3129  << tt_end-tt_start << std::endl;
3130  oomph_info.stream_pt()->flush();
3131  tt_start = TimingHelpers::timer();
3132  }
3133 
3134  //MemoryUsage::doc_memory_usage("after setup loops");
3135 
3136  // At this point we need to "synchronise" the nodes on halo(ed) elements
3137  // so that the processors_associated_with_data agrees for the same node
3138  // on all processors. Strategy: All local nodes have just had their
3139  // association recorded. Now loop over all haloed elements
3140  // and send the association of their nodes to the corresponding
3141  // halo processors where they update/augment the association of the
3142  // nodes of the corresponding halo elements.
3143 
3144  // Loop over all domains
3145  for (int d=0;d<n_proc;d++)
3146  {
3147  // Prepare vector to send/receive
3148  Vector<unsigned> processors_associated_with_data_on_other_proc;
3149 
3150  if (d!=my_rank)
3151  {
3152  // Communicate the processors associated with data on haloed elements
3153 
3154  // Get haloed elements
3155  Vector<GeneralisedElement*> haloed_elem_pt(this->haloed_element_pt(d));
3156 
3157  // Initialise counter for this haloed layer
3158  unsigned count_data=0;
3159 
3160  // Loop over haloed elements
3161  unsigned n_haloed_elem=haloed_elem_pt.size();
3162  for (unsigned e=0;e<n_haloed_elem;e++)
3163  {
3164  //Only nodes in finite elements
3165  FiniteElement* haloed_el_pt =
3166  dynamic_cast<FiniteElement*>(haloed_elem_pt[e]);
3167  if(haloed_el_pt!=0)
3168  {
3169  // Loop over nodes
3170  unsigned n_node=haloed_el_pt->nnode();
3171  for (unsigned j=0;j<n_node;j++)
3172  {
3173  Node* nod_pt=haloed_el_pt->node_pt(j);
3174 
3175  // Number of processors associated with this node
3176  unsigned n_assoc=processors_associated_with_data[nod_pt].size();
3177 
3178  // This number needs to be sent
3179  processors_associated_with_data_on_other_proc.push_back(n_assoc);
3180  count_data++;
3181 
3182  // Now add the process IDs associated to the vector to be sent
3183  std::set<unsigned> procs_set=
3184  processors_associated_with_data[nod_pt];
3185  for (std::set<unsigned>::iterator it=procs_set.begin();
3186  it!=procs_set.end();it++)
3187  {
3188  processors_associated_with_data_on_other_proc.push_back(*it);
3189  count_data++;
3190  }
3191  }
3192  }
3193  }
3194 
3195 
3196 
3197  // Send the information
3198  MPI_Send(&count_data,1,MPI_UNSIGNED,d,0,Comm_pt->mpi_comm());
3199  if (count_data!=0)
3200  {
3201  MPI_Send(&processors_associated_with_data_on_other_proc[0],count_data,
3202  MPI_UNSIGNED,d,1,Comm_pt->mpi_comm());
3203  }
3204  }
3205  else
3206  {
3207  // Receive the processors associated with data onto halo elements
3208  for (int dd=0;dd<n_proc;dd++)
3209  {
3210  if (dd!=my_rank) // (my_rank=d)
3211  {
3212  // We will be looping over the halo elements with process dd
3213  Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(dd));
3214  unsigned n_halo_elem=halo_elem_pt.size();
3215  unsigned count_data=0;
3216  MPI_Recv(&count_data,1,MPI_UNSIGNED,dd,0,Comm_pt->mpi_comm(),&status);
3217  if (count_data!=0)
3218  {
3219  processors_associated_with_data_on_other_proc.resize(count_data);
3220  MPI_Recv(&processors_associated_with_data_on_other_proc[0],
3221  count_data,MPI_UNSIGNED,dd,1,Comm_pt->mpi_comm(),&status);
3222 
3223  // Reset counter and loop through nodes on halo elements
3224  count_data=0;
3225  for (unsigned e=0;e<n_halo_elem;e++)
3226  {
3227  FiniteElement* halo_el_pt=
3228  dynamic_cast<FiniteElement*>(halo_elem_pt[e]);
3229  if(halo_el_pt!=0)
3230  {
3231  unsigned n_node=halo_el_pt->nnode();
3232  for (unsigned j=0;j<n_node;j++)
3233  {
3234  Node* nod_pt=halo_el_pt->node_pt(j);
3235 
3236  // Get number of processors associated with data that was sent
3237  unsigned n_assoc=
3238  processors_associated_with_data_on_other_proc[count_data];
3239  count_data++;
3240 
3241  for (unsigned i_assoc=0;i_assoc<n_assoc;i_assoc++)
3242  {
3243  // Get the process ID
3244  unsigned sent_domain=
3245  processors_associated_with_data_on_other_proc[count_data];
3246  count_data++;
3247 
3248  // Add it to this processor's list of IDs
3249  processors_associated_with_data[nod_pt].insert(sent_domain);
3250 
3251  // If the node is solid then add the ID to the solid data
3252  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
3253  if (solid_nod_pt!=0)
3254  {
3255  processors_associated_with_data
3256  [solid_nod_pt->variable_position_pt()].
3257  insert(sent_domain);
3258  }
3259  }
3260  }
3261  }
3262  }
3263  }
3264  }
3265  }
3266  }
3267  }
3268 
3270  {
3271  tt_end = TimingHelpers::timer();
3272  oomph_info
3273  << "Time for pt2pt send/recv in Mesh::classify_halo_and_haloed_nodes: "
3274  << tt_end-tt_start << std::endl;
3275  oomph_info.stream_pt()->flush();
3276  tt_start = TimingHelpers::timer();
3277  }
3278 
3279 
3280  //MemoryUsage::doc_memory_usage("after pt2pt send/recv");
3281 
3282  // Loop over all nodes on the present processor and put the highest-numbered
3283  // processor associated with each node "in charge" of the node
3284  unsigned nnod=this->nnode();
3285  for (unsigned j=0;j<nnod;j++)
3286  {
3287  Node* nod_pt=this->node_pt(j);
3288 
3289  // Reset halo status of node to false
3290  nod_pt->set_nonhalo();
3291 
3292  // If it's a SolidNode then the halo status of the data
3293  // associated with that must also be reset to false
3294  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
3295  if (solid_nod_pt!=0)
3296  {
3297  solid_nod_pt->variable_position_pt()->set_nonhalo();
3298  }
3299 
3300  // Now put the highest-numbered one in charge
3301  unsigned proc_max=0;
3302  std::set<unsigned> procs_set=processors_associated_with_data[nod_pt];
3303  for (std::set<unsigned>::iterator it=procs_set.begin();
3304  it!=procs_set.end();it++)
3305  {
3306  if (*it>proc_max) proc_max=*it;
3307  }
3308  processor_in_charge[nod_pt]=proc_max;
3309 
3310  // Do the same if we have a SolidNode
3311  if (solid_nod_pt!=0)
3312  {
3313  // Now put the highest-numbered one in charge
3314  unsigned proc_max_solid=0;
3315  std::set<unsigned> procs_set_solid=processors_associated_with_data
3316  [solid_nod_pt->variable_position_pt()];
3317  for (std::set<unsigned>::iterator it=procs_set_solid.begin();
3318  it!=procs_set_solid.end();it++)
3319  {
3320  if (*it>proc_max_solid) proc_max_solid=*it;
3321  }
3322  processor_in_charge[solid_nod_pt->variable_position_pt()]=proc_max_solid;
3323  }
3324  }
3325 
3326 
3327  // First stab at determining halo nodes. They are located on the halo
3328  // elements and the processor in charge differs from the
3329  // current processor
3330 
3331  // Only count nodes once (map is initialised to 0 = false)
3332  std::map<Node*,bool> done;
3333 
3334  // Loop over all processors
3335  for (int domain=0;domain<n_proc;domain++)
3336  {
3337  // Get vector of halo elements by copy operation
3338  Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(domain));
3339 
3340  // Loop over halo elements associated with this adjacent domain
3341  unsigned nelem=halo_elem_pt.size();
3342 
3343  for (unsigned e=0;e<nelem;e++)
3344  {
3345  // Get element
3346  GeneralisedElement* el_pt=halo_elem_pt[e];
3347 
3348  //Can if be cast to a finite element
3349  FiniteElement* finite_el_pt=dynamic_cast<FiniteElement*>(el_pt);
3350  if(finite_el_pt!=0)
3351  {
3352  //Loop over nodes
3353  unsigned nnod=finite_el_pt->nnode();
3354  for (unsigned j=0;j<nnod;j++)
3355  {
3356  Node* nod_pt=finite_el_pt->node_pt(j);
3357 
3358  // Have we done this node already?
3359  if (!done[nod_pt])
3360  {
3361  // Is the other processor/domain in charge of this node?
3362  int proc_in_charge=processor_in_charge[nod_pt];
3363 
3364  if (proc_in_charge!=my_rank)
3365  {
3366 
3367  // To keep the order of the nodes consistent with that
3368  // in the haloed node lookup scheme, only
3369  // allow it to be added when the current domain is in charge
3370  if (proc_in_charge==int(domain))
3371  {
3372  // Add it as being halo node whose non-halo counterpart
3373  // is located on processor proc_in_charge
3374  this->add_halo_node_pt(proc_in_charge,nod_pt);
3375 
3376  // The node itself needs to know it is a halo
3377  nod_pt->set_halo(proc_in_charge);
3378 
3379  // If it's a SolidNode then the data associated with that
3380  // must also be halo
3381  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
3382  if (solid_nod_pt!=0)
3383  {
3384  solid_nod_pt->variable_position_pt()->
3385  set_halo(proc_in_charge);
3386  }
3387 
3388  // We're done with this node
3389  done[nod_pt]=true;
3390  }
3391  }
3392  }
3393  }
3394 
3395  } //End of finite element case
3396 
3397  // Now make sure internal data on halo elements is also halo
3398  unsigned nintern_data = el_pt->ninternal_data();
3399  for (unsigned iintern=0;iintern<nintern_data;iintern++)
3400  {
3401  el_pt->internal_data_pt(iintern)->set_halo(domain);
3402  }
3403  }
3404  }
3405 
3406 
3407  // First stab at determining haloed nodes. They are located on the haloed
3408  // elements and the processor in charge is the current processor
3409 
3410  // Loop over processors that share haloes with the current one
3411  for (int domain=0;domain<n_proc;domain++)
3412  {
3413  // Only count nodes once (map is initialised to 0 = false)
3414  std::map<Node*,bool> node_done;
3415 
3416  // Get vector of haloed elements by copy operation
3417  Vector<GeneralisedElement*> haloed_elem_pt(this->haloed_element_pt(domain));
3418 
3419  // Loop over haloed elements associated with this adjacent domain
3420  unsigned nelem=haloed_elem_pt.size();
3421 
3422  for (unsigned e=0;e<nelem;e++)
3423  {
3424  // Get element
3425  GeneralisedElement* el_pt=haloed_elem_pt[e];
3426 
3427  //Can it be cast to a finite element
3428  FiniteElement* finite_el_pt=dynamic_cast<FiniteElement*>(el_pt);
3429  if(finite_el_pt!=0)
3430  {
3431  //Loop over nodes
3432  unsigned nnod=finite_el_pt->nnode();
3433  for (unsigned j=0;j<nnod;j++)
3434  {
3435  Node* nod_pt=finite_el_pt->node_pt(j);
3436 
3437  // Have we done this node already?
3438  if (!node_done[nod_pt])
3439  {
3440  // Is the current processor/domain in charge of this node?
3441  int proc_in_charge=processor_in_charge[nod_pt];
3442 
3443  if (proc_in_charge==my_rank)
3444  {
3445  // Add it as being haloed from specified domain
3446  this->add_haloed_node_pt(domain,nod_pt);
3447 
3448  // We're done with this node
3449  node_done[nod_pt]=true;
3450  }
3451  }
3452 
3453  }
3454  }
3455  }
3456  }
3457 
3458 
3460  {
3461  tt_end = TimingHelpers::timer();
3462  oomph_info
3463  << "Time for first classific in Mesh::classify_halo_and_haloed_nodes: "
3464  << tt_end-tt_start << std::endl;
3465  oomph_info.stream_pt()->flush();
3466  tt_start = TimingHelpers::timer();
3467  }
3468 
3469  //MemoryUsage::doc_memory_usage("after first classific");
3470 
3471 
3472  // Find any overlooked halo nodes: These are any nodes on the halo/haloed
3473  // elements (i.e. precisely the nodes currently contained in the shared
3474  // node scheme) that have not been classified as haloes (yet) though they
3475  // should have been because another processor is in charge of them.
3476  // This arises when the "overlooked halo node" is not part of the
3477  // halo/haloed element lookup scheme between the current processor
3478  // and the processor that holds the non-halo counterpart. E.g. we're
3479  // on proc 3. A node at the very edge of its halo layer also exists
3480  // on procs 0 and 1 with 1 being "in charge". However, the node in
3481  // question is not part of the halo/haloed element lookup scheme between
3482  // processor 1 and 3 so in the classification performed above, we never
3483  // visit it so it's overlooked. The code below rectifies this by going
3484  // through the intermediate processor (here proc 0) that contains the node in
3485  // lookup schemes with the halo processor (here proc 3, this one) and the one
3486  // that contains the non-halo counterpart (here proc 1).
3487 
3488 
3489  // Counter for number of overlooked halos (if there aren't any we don't
3490  // need any comms below)
3491  unsigned n_overlooked_halo=0;
3492 
3493  // Record previously overlooked halo nodes so they can be
3494  // added to the shared node lookup scheme (in a consistent order) below
3495  Vector<Vector<Node*> > over_looked_halo_node_pt(n_proc);
3496 
3497  // Record previously overlooked haloed nodes so they can be
3498  // added to the shared node lookup scheme (in a consistent order) below
3499  Vector<Vector<Node*> > over_looked_haloed_node_pt(n_proc);
3500 
3501  // Data to be sent to each processor
3502  Vector<int> send_n(n_proc,0);
3503 
3504  // Storage for all values to be sent to all processors
3505  Vector<int> send_data;
3506 
3507  // Start location within send_data for data to be sent to each processor
3508  Vector<int> send_displacement(n_proc,0);
3509 
3510  // Check missing ones
3511  for (int domain=0;domain<n_proc;domain++)
3512  {
3513  //Set the offset for the current processor
3514  send_displacement[domain] = send_data.size();
3515 
3516  //Don't bother to do anything if the processor in the loop is the
3517  //current processor
3518  if(domain!=my_rank)
3519  {
3520  unsigned nnod=nshared_node(domain);
3521  for (unsigned j=0;j<nnod;j++)
3522  {
3523  Node* nod_pt=shared_node_pt(domain,j);
3524 
3525  // Is a different-numbered processor/domain in charge of this node?
3526  int proc_in_charge=processor_in_charge[nod_pt];
3527  if ((proc_in_charge!=my_rank)&&!(nod_pt->is_halo()))
3528  {
3529  // Add it as being halo node whose non-halo counterpart
3530  // is located on processor proc_in_charge
3531  this->add_halo_node_pt(proc_in_charge,nod_pt);
3532 
3533  // We have another one...
3534  n_overlooked_halo++;
3535  over_looked_halo_node_pt[proc_in_charge].push_back(nod_pt);
3536 
3537  // The node itself needs to know it is a halo
3538  nod_pt->set_halo(proc_in_charge);
3539 
3540  // If it's a SolidNode then the data associated with that
3541  // must also be halo
3542  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
3543  if (solid_nod_pt!=0)
3544  {
3545  solid_nod_pt->variable_position_pt()->set_halo(proc_in_charge);
3546  }
3547 
3548  // Send shared node ID and processor in charge info to
3549  // "intermediate" processor (i.e. the processor that has
3550  // the lookup schemes to talk to both
3551  send_data.push_back(j);
3552  send_data.push_back(proc_in_charge);
3553  }
3554  }
3555  }
3556 
3557  // End of data
3558  send_data.push_back(-1);
3559 
3560  //Find the number of data added to the vector
3561  send_n[domain] = send_data.size() - send_displacement[domain];
3562  }
3563 
3564 
3565  // Check if any processor has stumbled across overlooked halos
3566  // (if not we can omit the comms below)
3567  unsigned global_max_n_overlooked_halo=0;
3568  MPI_Allreduce(&n_overlooked_halo,&global_max_n_overlooked_halo,1,
3569  MPI_UNSIGNED,MPI_MAX,
3570  Comm_pt->mpi_comm());
3571 
3572 
3573  oomph_info << "Global max number of overlooked haloes: "
3574  << global_max_n_overlooked_halo << std::endl;
3575 
3577  {
3578  tt_end = TimingHelpers::timer();
3579  oomph_info
3580  <<"Time for setup 1st alltoalls in Mesh::classify_halo_and_haloed_nodes: "
3581  << tt_end-tt_start << std::endl;
3582  oomph_info.stream_pt()->flush();
3583  tt_start = TimingHelpers::timer();
3584  }
3585 
3586  //MemoryUsage::doc_memory_usage("after setup 1st alltoalls");
3587 
3588  // Any comms needed?
3589  if (global_max_n_overlooked_halo>0)
3590  {
3591 
3592 
3593 
3594  //Storage for the number of data to be received from each processor
3595  Vector<int> receive_n(n_proc,0);
3596 
3597  //Now send numbers of data to be sent between all processors
3598  MPI_Alltoall(&send_n[0],1,MPI_INT,&receive_n[0],1,MPI_INT,
3599  Comm_pt->mpi_comm());
3600 
3601 
3602 
3604  {
3605  tt_end = TimingHelpers::timer();
3606  oomph_info
3607  << "Time for 1st alltoall in Mesh::classify_halo_and_haloed_nodes: "
3608  << tt_end-tt_start << std::endl;
3609  oomph_info.stream_pt()->flush();
3610  tt_start = TimingHelpers::timer();
3611  }
3612 
3613  //MemoryUsage::doc_memory_usage("after 1st alltoall");
3614 
3615 
3616  //We now prepare the data to be received
3617  //by working out the displacements from the received data
3618  Vector<int> receive_displacement(n_proc,0);
3619  int receive_data_count=0;
3620  for(int rank=0;rank<n_proc;++rank)
3621  {
3622  //Displacement is number of data received so far
3623  receive_displacement[rank] = receive_data_count;
3624  receive_data_count += receive_n[rank];
3625  }
3626 
3627  //Now resize the receive buffer for all data from all processors
3628  //Make sure that it has a size of at least one
3629  if(receive_data_count==0) {++receive_data_count;}
3630  Vector<int> receive_data(receive_data_count);
3631 
3632  //Make sure that the send buffer has size at least one
3633  //so that we don't get a segmentation fault
3634  if(send_data.size()==0) {send_data.resize(1);}
3635 
3636  //Now send the data between all the processors
3637  MPI_Alltoallv(&send_data[0],&send_n[0],&send_displacement[0],
3638  MPI_INT,
3639  &receive_data[0],&receive_n[0],
3640  &receive_displacement[0],
3641  MPI_INT,
3642  Comm_pt->mpi_comm());
3643 
3644 
3645 
3647  {
3648  tt_end = TimingHelpers::timer();
3649  oomph_info
3650  << "Time for 2nd alltoall in Mesh::classify_halo_and_haloed_nodes: "
3651  << tt_end-tt_start << std::endl;
3652  oomph_info.stream_pt()->flush();
3653  tt_start = TimingHelpers::timer();
3654  }
3655 
3656  //MemoryUsage::doc_memory_usage("after 2nd alltoall");
3657 
3658  // Provide storage for data to be sent to processor that used to be
3659  // in charge
3660  Vector<Vector<int> > send_data_for_proc_in_charge(n_proc);
3661 
3662  //Now use the received data
3663  for (int send_rank=0;send_rank<n_proc;send_rank++)
3664  {
3665  //Don't bother to do anything for the processor corresponding to the
3666  //current processor or if no data were received from this processor
3667  if((send_rank != my_rank) && (receive_n[send_rank] != 0))
3668  {
3669  //Counter for the data within the large array
3670  unsigned count=receive_displacement[send_rank];
3671 
3672  // Unpack until we reach "end of data" indicator (-1)
3673  while(true)
3674  {
3675 
3676  //Read next entry
3677  int next_one=receive_data[count++];
3678 
3679  if (next_one==-1)
3680  {
3681  break;
3682  }
3683  else
3684  {
3685  // Shared halo node number in lookup scheme between intermediate
3686  // (i.e. this) processor and the one that has the overlooked halo
3687  unsigned j=unsigned(next_one);
3688 
3689  // Processor in charge:
3690  unsigned proc_in_charge=unsigned(receive_data[count++]);
3691 
3692  // Find actual node from shared node lookup scheme
3693  Node* nod_pt=shared_node_pt(send_rank,j);
3694 
3695 
3696  // Note: This search seems relatively cheap
3697  // and in the tests done, did not benefit
3698  // from conversion to map-based search
3699  // as in
3700  // TreeBasedRefineableMeshBase::synchronise_hanging_nodes()
3701 
3702 
3703  // Locate its index in lookup scheme with proc in charge
3704  bool found=false;
3705  unsigned nnod=nshared_node(proc_in_charge);
3706  for (unsigned jj=0;jj<nnod;jj++)
3707  {
3708  if (nod_pt==shared_node_pt(proc_in_charge,jj))
3709  {
3710  found=true;
3711 
3712  // Shared node ID in lookup scheme with intermediate (i.e. this)
3713  // processor
3714  send_data_for_proc_in_charge[proc_in_charge].push_back(jj);
3715 
3716  // Processor that holds the overlooked halo node
3717  send_data_for_proc_in_charge[proc_in_charge].push_back(
3718  send_rank);
3719 
3720  break;
3721  }
3722  }
3723  if (!found)
3724  {
3725 
3726  std::ostringstream error_stream;
3727  error_stream << "Failed to find node that is shared node " << j
3728  << " (with processor " << send_rank
3729  << ") \n in shared node lookup scheme with processor "
3730  << proc_in_charge << " which is in charge.\n";
3731  throw OomphLibError(error_stream.str(),
3732  OOMPH_CURRENT_FUNCTION,
3733  OOMPH_EXCEPTION_LOCATION);
3734  }
3735  }
3736  }
3737 
3738  }
3739  } //End of data is received
3740 
3741 
3743  {
3744  tt_end = TimingHelpers::timer();
3745  oomph_info
3746  << "Time for 1st setup 3rd alltoall in Mesh::classify_halo_and_haloed_nodes: "
3747  << tt_end-tt_start << std::endl;
3748  oomph_info.stream_pt()->flush();
3749  tt_start = TimingHelpers::timer();
3750  }
3751 
3752  //MemoryUsage::doc_memory_usage("after 1st setup for 3rd alltoall");
3753 
3754  // Data to be sent to each processor
3755  Vector<int> all_send_n(n_proc,0);
3756 
3757  // Storage for all values to be sent to all processors
3758  Vector<int> all_send_data;
3759 
3760  // Start location within send_data for data to be sent to each processor
3761  Vector<int> all_send_displacement(n_proc,0);
3762 
3763  // Collate data
3764  for (int domain=0;domain<n_proc;domain++)
3765  {
3766  //Set the offset for the current processor
3767  all_send_displacement[domain] = all_send_data.size();
3768 
3769  //Don't bother to do anything if the processor in the loop is the
3770  //current processor
3771  if(domain!=my_rank)
3772  {
3773  unsigned n=send_data_for_proc_in_charge[domain].size();
3774  for (unsigned j=0;j<n;j++)
3775  {
3776  all_send_data.push_back(
3777  send_data_for_proc_in_charge[domain][j]);
3778  }
3779  }
3780 
3781  // End of data
3782  all_send_data.push_back(-1);
3783 
3784  //Find the number of data added to the vector
3785  all_send_n[domain]=all_send_data.size()-all_send_displacement[domain];
3786  }
3787 
3788 
3790  {
3791  tt_end = TimingHelpers::timer();
3792  oomph_info
3793  << "Time for 2nd setup 3rd alltoall in Mesh::classify_halo_and_haloed_nodes: "
3794  << tt_end-tt_start << std::endl;
3795  oomph_info.stream_pt()->flush();
3796  tt_start = TimingHelpers::timer();
3797  }
3798 
3799  //MemoryUsage::doc_memory_usage("after 2nd setup 3rd alltoall");
3800 
3801  //Storage for the number of data to be received from each processor
3802  Vector<int> all_receive_n(n_proc,0);
3803 
3804  //Now send numbers of data to be sent between all processors
3805  MPI_Alltoall(&all_send_n[0],1,MPI_INT,&all_receive_n[0],1,MPI_INT,
3806  Comm_pt->mpi_comm());
3807 
3808 
3809 
3811  {
3812  tt_end = TimingHelpers::timer();
3813  oomph_info
3814  << "Time for 3rd alltoall in Mesh::classify_halo_and_haloed_nodes: "
3815  << tt_end-tt_start << std::endl;
3816  oomph_info.stream_pt()->flush();
3817  tt_start = TimingHelpers::timer();
3818  }
3819 
3820  //MemoryUsage::doc_memory_usage("after 3rd alltoall");
3821 
3822  //We now prepare the data to be received
3823  //by working out the displacements from the received data
3824  Vector<int> all_receive_displacement(n_proc,0);
3825  int all_receive_data_count=0;
3826 
3827  for(int rank=0;rank<n_proc;++rank)
3828  {
3829  //Displacement is number of data received so far
3830  all_receive_displacement[rank] = all_receive_data_count;
3831  all_receive_data_count += all_receive_n[rank];
3832  }
3833 
3834  //Now resize the receive buffer for all data from all processors
3835  //Make sure that it has a size of at least one
3836  if (all_receive_data_count==0) {++all_receive_data_count;}
3837  Vector<int> all_receive_data(all_receive_data_count);
3838 
3839  //Make sure that the send buffer has size at least one
3840  //so that we don't get a segmentation fault
3841  if(all_send_data.size()==0) {all_send_data.resize(1);}
3842 
3843  //Now send the data between all the processors
3844  MPI_Alltoallv(&all_send_data[0],&all_send_n[0],&all_send_displacement[0],
3845  MPI_INT,
3846  &all_receive_data[0],&all_receive_n[0],
3847  &all_receive_displacement[0],
3848  MPI_INT,
3849  Comm_pt->mpi_comm());
3850 
3851 
3853  {
3854  tt_end = TimingHelpers::timer();
3855  oomph_info
3856  << "Time for 4th alltoall in Mesh::classify_halo_and_haloed_nodes: "
3857  << tt_end-tt_start << std::endl;
3858  oomph_info.stream_pt()->flush();
3859  tt_start = TimingHelpers::timer();
3860  }
3861 
3862  //MemoryUsage::doc_memory_usage("after 4th alltoall");
3863 
3864 
3865  //Now use the received data
3866  for (int send_rank=0;send_rank<n_proc;send_rank++)
3867  {
3868  //Don't bother to do anything for the processor corresponding to the
3869  //current processor or if no data were received from this processor
3870  if((send_rank != my_rank) && (all_receive_n[send_rank] != 0))
3871  {
3872  //Counter for the data within the large array
3873  unsigned count=all_receive_displacement[send_rank];
3874 
3875  // Unpack until we reach "end of data" indicator (-1)
3876  while(true)
3877  {
3878 
3879  //Read next entry
3880  int next_one=all_receive_data[count++];
3881 
3882  if (next_one==-1)
3883  {
3884  break;
3885  }
3886  else
3887  {
3888  // Shared node ID in lookup scheme with intermediate (sending)
3889  // processor
3890  unsigned j=unsigned(next_one);
3891 
3892  // Get pointer to previously overlooked halo
3893  Node* nod_pt=shared_node_pt(send_rank,j);
3894 
3895  // Proc where overlooked halo is
3896  unsigned proc_with_overlooked_halo=all_receive_data[count++];
3897 
3898  // Add it as being haloed from specified domain
3899  this->add_haloed_node_pt(proc_with_overlooked_halo,nod_pt);
3900 
3901  // Record new haloed node so it an be added to the shared
3902  // node lookup scheme (in a consistent order) below.
3903  over_looked_haloed_node_pt[proc_with_overlooked_halo].
3904  push_back(nod_pt);
3905  }
3906  }
3907  }
3908  } //End of data is received
3909 
3911  {
3912  tt_end = TimingHelpers::timer();
3913  oomph_info
3914  << "Time for postproc 4th alltoall in Mesh::classify_halo_and_haloed_nodes: "
3915  << tt_end-tt_start << std::endl;
3916  oomph_info.stream_pt()->flush();
3917  tt_start = TimingHelpers::timer();
3918  }
3919 
3920  //MemoryUsage::doc_memory_usage("after postprocess 4th alltoall");
3921 
3922  // Now add previously overlooked halo/haloed nodes to shared node
3923  // lookup scheme in consistent order
3924  for (int d=0;d<n_proc;d++)
3925  {
3926  // For all domains lower than the current domain: Do halos first
3927  // then haloed, to ensure correct order in lookup scheme from
3928  // the other side
3929  if (d<my_rank)
3930  {
3931  unsigned nnod=over_looked_halo_node_pt[d].size();
3932  for (unsigned j=0;j<nnod;j++)
3933  {
3934  this->add_shared_node_pt(d,over_looked_halo_node_pt[d][j]);
3935  }
3936  nnod=over_looked_haloed_node_pt[d].size();
3937  for (unsigned j=0;j<nnod;j++)
3938  {
3939  this->add_shared_node_pt(d,over_looked_haloed_node_pt[d][j]);
3940  }
3941  }
3942  else if (d>my_rank)
3943  {
3944 
3945  unsigned nnod=over_looked_haloed_node_pt[d].size();
3946  for (unsigned j=0;j<nnod;j++)
3947  {
3948  this->add_shared_node_pt(d,over_looked_haloed_node_pt[d][j]);
3949  }
3950  nnod=over_looked_halo_node_pt[d].size();
3951  for (unsigned j=0;j<nnod;j++)
3952  {
3953  this->add_shared_node_pt(d,over_looked_halo_node_pt[d][j]);
3954  }
3955  }
3956  }
3957 
3958 
3959  // Doc stats
3960  if (report_stats)
3961  {
3962  // Report total number of halo(ed) and shared nodes for this process
3963  oomph_info << "BEFORE SYNCHRONISE SHARED NODES Processor " << my_rank
3964  << " holds " << this->nnode()
3965  << " nodes of which " << this->nhalo_node()
3966  << " are halo nodes \n while " << this->nhaloed_node()
3967  << " are haloed nodes, and " << this->nshared_node()
3968  << " are shared nodes." << std::endl;
3969 
3970  // Report number of halo(ed) and shared nodes with each domain
3971  // from the current process
3972  for (int iproc=0;iproc<n_proc;iproc++)
3973  {
3974  // Get vector of halo elements by copy operation
3975  Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(iproc));
3976  Vector<GeneralisedElement*> haloed_elem_pt(
3977  this->haloed_element_pt(iproc));
3978  oomph_info << "With process " << iproc << ", there are "
3979  << this->nhalo_node(iproc) << " halo nodes, and " << std::endl
3980  << this->nhaloed_node(iproc) << " haloed nodes, and "
3981  << this->nshared_node(iproc) << " shared nodes" << std::endl
3982  << halo_elem_pt.size() << " halo elements and "
3983  << haloed_elem_pt.size() << " haloed elements\n";
3984  }
3985  }
3986 
3987 // // Doc stats
3988 // if (report_stats)
3989 // {
3990 // // Report total number of halo(ed) and shared nodes for this process
3991 // oomph_info << "BEFORE SYNCHRONISE SHARED NODES Processor " << my_rank
3992 // << " holds " << this->nnode()
3993 // << " nodes of which " << this->nhalo_node()
3994 // << " are halo nodes \n while " << this->nhaloed_node()
3995 // << " are haloed nodes, and " << this->nshared_node()
3996 // << " are shared nodes." << std::endl;
3997 
3998 // // Report number of halo(ed) and shared nodes with each domain
3999 // // from the current process
4000 // for (int iproc=0;iproc<n_proc;iproc++)
4001 // {
4002 // // Get vector of halo elements by copy operation
4003 // Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(iproc));
4004 // Vector<GeneralisedElement*> haloed_elem_pt(
4005 // this->haloed_element_pt(iproc));
4006 // oomph_info << "With process " << iproc << ", there are "
4007 // << this->nhalo_node(iproc) << " halo nodes, and " << std::endl
4008 // << this->nhaloed_node(iproc) << " haloed nodes, and "
4009 // << this->nshared_node(iproc) << " shared nodes" << std::endl
4010 // << halo_elem_pt.size() << " halo elements and "
4011 // << haloed_elem_pt.size() << " haloed elements\n";
4012 // }
4013 // }
4014 
4015  } // end if comms reqd because we encountered overlooked halo elements
4016 
4017 
4018  //MemoryUsage::doc_memory_usage("before sync halo nodes");
4019 
4020  // Synchronise shared nodes
4021  synchronise_shared_nodes(report_stats);
4022 
4023  //MemoryUsage::doc_memory_usage("after sync halo nodes");
4024 
4025 #ifdef PARANOID
4026  // Has some twit filled in haloed nodes with own process?!
4027  if (Haloed_node_pt[my_rank].size()!=0)
4028  {
4029  throw OomphLibError(
4030  "Processor has haloed nodes with itself! Something's gone wrong!",
4031  OOMPH_CURRENT_FUNCTION,
4032  OOMPH_EXCEPTION_LOCATION);
4033  }
4034  // Has some twit filled in halo nodes with own process?!
4035  if (Halo_node_pt[my_rank].size()!=0)
4036  {
4037  throw OomphLibError(
4038  "Processor has halo nodes with itself! Something's gone wrong!",
4039  OOMPH_CURRENT_FUNCTION,
4040  OOMPH_EXCEPTION_LOCATION);
4041  }
4042  // Has some twit filled in root haloed elements with own process?!
4043  if (Root_haloed_element_pt[my_rank].size()!=0)
4044  {
4045  throw OomphLibError(
4046  "Processor has root haloed elements with itself! Something's gone wrong!",
4047  OOMPH_CURRENT_FUNCTION,
4048  OOMPH_EXCEPTION_LOCATION);
4049  }
4050  // Has some twit filled in root halo elements with own process?!
4051  if (Root_halo_element_pt[my_rank].size()!=0)
4052  {
4053  throw OomphLibError(
4054  "Processor has root halo elements with itself! Something's gone wrong!",
4055  OOMPH_CURRENT_FUNCTION,
4056  OOMPH_EXCEPTION_LOCATION);
4057  }
4058 #endif
4059 
4060  // Doc stats
4061  if (report_stats)
4062  {
4063  // Report total number of halo(ed) and shared nodes for this process
4064  oomph_info << "Processor " << my_rank
4065  << " holds " << this->nnode()
4066  << " nodes of which " << this->nhalo_node()
4067  << " are halo nodes \n while " << this->nhaloed_node()
4068  << " are haloed nodes, and " << this->nshared_node()
4069  << " are shared nodes." << std::endl;
4070 
4071  // Report number of halo(ed) and shared nodes with each domain
4072  // from the current process
4073  for (int iproc=0;iproc<n_proc;iproc++)
4074  {
4075  // Get vector of halo elements by copy operation
4076  Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(iproc));
4077  Vector<GeneralisedElement*> haloed_elem_pt(
4078  this->haloed_element_pt(iproc));
4079  oomph_info << "With process " << iproc << ", there are "
4080  << this->nhalo_node(iproc) << " halo nodes, and " << std::endl
4081  << this->nhaloed_node(iproc) << " haloed nodes, and "
4082  << this->nshared_node(iproc) << " shared nodes" << std::endl
4083  << halo_elem_pt.size() << " halo elements and "
4084  << haloed_elem_pt.size() << " haloed elements\n";
4085  }
4086  }
4087 
4088 
4089  //MemoryUsage::doc_memory_usage("before resize halo nodes");
4090 
4091  // Now resize halo nodes if required (can be over-ruled from the outside
4092  // either by user (for (risky!) efficienty saving) or from overloaded
4093  // version of classify_... in refineable version of that function
4094  // where resize_halo_nodes() is called after synchronising hanging nodes.
4096  {
4098  }
4099 
4100  //MemoryUsage::doc_memory_usage("after resize halo nodes");
4101 
4103  {
4104  double t_end = TimingHelpers::timer();
4105  oomph_info << "Total time for Mesh::classify_halo_and_halo_nodes(): "
4106  << t_end-t_start << std::endl;
4107  oomph_info.stream_pt()->flush();
4108  }
4109 
4110 // MemoryUsage::doc_memory_usage(
4111 // "at end of Mesh::classify_halo_and_halo_nodes()");
4112 
4113 }
4114 
4115 
4116 //========================================================================
4117 /// Helper function that resizes halo nodes to the same
4118 /// size as their non-halo counterparts if required. (A discrepancy
4119 /// can arise if a FaceElement that introduces additional unknowns
4120 /// are attached to a bulk element that shares a node with a haloed element.
4121 /// In that case the joint node between haloed and non-haloed element
4122 /// is resized on that processor but not on the one that holds the
4123 /// halo counterpart (because no FaceElement is attached to the halo
4124 /// element)
4125 //=========================================================================
4127 {
4128 
4129 
4130  double t_start = 0.0;
4132  {
4133  t_start=TimingHelpers::timer();
4134  }
4135 
4136  MPI_Status status;
4137 
4138  // Nuffink needs to be done if mesh isn't distributed
4139  if (is_mesh_distributed())
4140  {
4141 
4142  // Storage for current processor and number of processors
4143  int n_proc=Comm_pt->nproc();
4144  int my_rank=Comm_pt->my_rank();
4145 
4146  // Loop over domains on which non-halo counter parts of my halo nodes live
4147  for (int d=0;d<n_proc;d++)
4148  {
4149  // On current processor: Receive data for my halo nodes with proc d.
4150  // Elsewhere: Send haloed data with proc d.
4151  if (d==my_rank)
4152  {
4153  // Loop over domains that hold non-halo counterparts of my halo nodes
4154  for (int dd=0;dd<n_proc;dd++)
4155  {
4156  // Don't talk to yourself
4157  if (dd!=d)
4158  {
4159  // How many of my nodes are halo nodes whose non-halo
4160  // counterpart is located on processor dd?
4161  int nnod_halo=this->nhalo_node(dd);
4162  int nnod_ext_halo=this->nexternal_halo_node(dd);
4163  if ((nnod_halo+nnod_ext_halo)!=0)
4164  {
4165  // Receive from processor dd number of haloed nodes (entry 0),
4166  // external haloed nodes (entry 1) and total number of
4167  // unsigneds to be sent below (entry 2)
4168  Vector<int> tmp(3);
4169  MPI_Recv(&tmp[0],3,MPI_INT,dd,0,Comm_pt->mpi_comm(),&status);
4170 
4171 #ifdef PARANOID
4172  // Check that number of halo/haloed nodes match
4173  int nnod_haloed=tmp[0];
4174  if (nnod_haloed!=nnod_halo)
4175  {
4176  std::ostringstream error_message;
4177  error_message
4178  << "Clash in numbers of halo and haloed nodes "
4179  << std::endl;
4180  error_message
4181  << " between procs "
4182  << dd << " and " << d << ": "
4183  << nnod_haloed << " " << nnod_halo << std::endl;
4184  throw OomphLibError(error_message.str(),
4185  OOMPH_CURRENT_FUNCTION,
4186  OOMPH_EXCEPTION_LOCATION);
4187  }
4188 
4189  // Check that number of external halo/haloed nodes match
4190  int nnod_ext_haloed=tmp[1];
4191  if (nnod_ext_haloed!=nnod_ext_halo)
4192  {
4193  std::ostringstream error_message;
4194  error_message
4195  << "Clash in numbers of external halo and haloed nodes "
4196  << std::endl;
4197  error_message
4198  << " between procs "
4199  << dd << " and " << d << ": "
4200  << nnod_ext_haloed << " " << nnod_ext_halo << std::endl;
4201  throw OomphLibError(error_message.str(),
4202  OOMPH_CURRENT_FUNCTION,
4203  OOMPH_EXCEPTION_LOCATION);
4204  }
4205 #endif
4206 
4207  // How many unsigneds are we about to receive
4208  unsigned n_rec=tmp[2];
4209 
4210  // Get strung-together data from other proc
4211  Vector<unsigned> unsigned_rec_data(n_rec);
4212  MPI_Recv(&unsigned_rec_data[0],n_rec,MPI_UNSIGNED,dd,
4213  0,Comm_pt->mpi_comm(),&status);
4214 
4215  // Step through the flat-packed unsigneds
4216  unsigned count=0;
4217 
4218  // Normal and external halo nodes
4219  for (unsigned loop=0;loop<2;loop++)
4220  {
4221  unsigned hi_nod=nnod_halo;
4222  if (loop==1)
4223  {
4224  hi_nod=nnod_ext_halo;
4225  }
4226  for (int j=0;j<int(hi_nod);j++)
4227  {
4228  // Which node are we dealing with
4229  Node* nod_pt=0;
4230  if (loop==0)
4231  {
4232  nod_pt=this->halo_node_pt(dd,j);
4233  }
4234  else
4235  {
4236  nod_pt=this->external_halo_node_pt(dd,j);
4237  }
4238 
4239  // How many values do we have locally?
4240  unsigned nval_local=nod_pt->nvalue();
4241 
4242  // Read number of values on other side
4243  unsigned nval_other=unsigned_rec_data[count++];
4244 
4245  if (nval_local!=nval_other)
4246  {
4247  nod_pt->resize(nval_other);
4248  }
4249 
4250  // Read number of entries in resize map
4251  unsigned nentry=unsigned_rec_data[count++];
4252  if (nentry!=0)
4253  {
4254  // Is current node a boundary node?
4255  BoundaryNodeBase* bnod_pt=
4256  dynamic_cast<BoundaryNodeBase*>(nod_pt);
4257 #ifdef PARANOID
4258  if (bnod_pt==0)
4259  {
4260  throw OomphLibError(
4261  "Failed to cast node to boundary node even though we've received data for boundary node",
4262  OOMPH_CURRENT_FUNCTION,
4263  OOMPH_EXCEPTION_LOCATION);
4264  }
4265 #endif
4266 
4267  // Create storage for map if it doesn't already exist
4268  bool already_existed=true;
4269  if(bnod_pt->
4270  index_of_first_value_assigned_by_face_element_pt()==0)
4271  {
4272  bnod_pt->
4273  index_of_first_value_assigned_by_face_element_pt()=
4274  new std::map<unsigned, unsigned>;
4275  already_existed=false;
4276  }
4277 
4278  // Get pointer to the map of indices associated with
4279  // additional values created by face elements
4280  std::map<unsigned, unsigned>* map_pt=
4282 
4283  // Loop over number of entries in map (as received)
4284  for (unsigned i=0;i<nentry;i++)
4285  {
4286  // Read out pairs...
4287  unsigned first_received=unsigned_rec_data[count++];
4288  unsigned second_received=unsigned_rec_data[count++];
4289 
4290  // If it exists check that values are consistent:
4291  if (already_existed)
4292  {
4293 #ifdef PARANOID
4294  if ((*map_pt)[first_received]!=second_received)
4295  {
4296  std::ostringstream error_message;
4297  error_message
4298  << "Existing map entry for map entry "
4299  << i << " for node located at ";
4300  unsigned n=nod_pt->ndim();
4301  for (unsigned ii=0;ii<n;ii++)
4302  {
4303  error_message << nod_pt->position(ii) << " ";
4304  }
4305  error_message
4306  << "Key: " << first_received << " "
4307  << "Local value: " << (*map_pt)[first_received] << " "
4308  << "Received value: " << second_received << std::endl;
4309  throw OomphLibError(error_message.str(),
4310  OOMPH_CURRENT_FUNCTION,
4311  OOMPH_EXCEPTION_LOCATION);
4312  }
4313 #endif
4314  }
4315  // Else assign
4316  else
4317  {
4318  (*map_pt)[first_received]=second_received;
4319  }
4320  }
4321  }
4322  }
4323  }
4324  }
4325  }
4326  }
4327  }
4328  // Send my haloed nodes whose halo counterparts are located on processor d
4329  else
4330  {
4331  // Storage for number of haloed nodes (entry 0), number of external
4332  // haloed nodes (entry 1) and total number of unsigneds to be sent below
4333  // (entry 2)
4334  Vector<int> tmp(3);
4335  int nnod_haloed=this->nhaloed_node(d);
4336  tmp[0]=nnod_haloed;
4337  int nnod_ext_haloed=this->nexternal_haloed_node(d);
4338  tmp[1]=nnod_ext_haloed;
4339  if ((nnod_haloed+nnod_ext_haloed)!=0)
4340  {
4341  // Now string together the data
4342  Vector<unsigned> unsigned_send_data;
4343 
4344  // Normal and external haloed nodes
4345  for (unsigned loop=0;loop<2;loop++)
4346  {
4347  unsigned hi_nod=nnod_haloed;
4348  if (loop==1)
4349  {
4350  hi_nod=nnod_ext_haloed;
4351  }
4352  for (int j=0;j<int(hi_nod);j++)
4353  {
4354  // Which node are we dealing with?
4355  Node* nod_pt=0;
4356  if (loop==0)
4357  {
4358  nod_pt=this->haloed_node_pt(d,j);
4359  }
4360  else
4361  {
4362  nod_pt=this->external_haloed_node_pt(d,j);
4363  }
4364 
4365  // Add number of values of node
4366  unsigned_send_data.push_back(nod_pt->nvalue());
4367 
4368  // Get pointer to the map of indices associated with
4369  // additional values created by face elements
4370 
4371  // Is it a boundary node?
4372  BoundaryNodeBase* bnod_pt=dynamic_cast<BoundaryNodeBase*>(nod_pt);
4373  if (bnod_pt==0)
4374  {
4375  // Not a boundary node -- there are zero map entries to follow
4376  unsigned_send_data.push_back(0);
4377  }
4378  else
4379  {
4380  // It's a boundary node: Check if it's been resized
4381  std::map<unsigned, unsigned>* map_pt=
4383 
4384  // No additional values created -- there are zero map entries
4385  // to follow
4386  if (map_pt==0)
4387  {
4388  unsigned_send_data.push_back(0);
4389  }
4390  // Created additional values
4391  else
4392  {
4393  // How many map entries were there
4394  unsigned_send_data.push_back(map_pt->size());
4395 
4396  // Loop over entries in map and add to send data
4397  for (std::map<unsigned, unsigned>::iterator p=
4398  map_pt->begin();
4399  p!=map_pt->end();p++)
4400  {
4401  unsigned_send_data.push_back((*p).first);
4402  unsigned_send_data.push_back((*p).second);
4403  }
4404  }
4405  }
4406  }
4407  }
4408 
4409  // How many values are there in total?
4410  int n_send=unsigned_send_data.size();
4411  tmp[2]=n_send;
4412 
4413  // Send the counts across to the other processor
4414  MPI_Send(&tmp[0],3,MPI_INT,d,0,Comm_pt->mpi_comm());
4415 
4416  // Send it across to the processor
4417  MPI_Send(&unsigned_send_data[0],n_send,MPI_UNSIGNED,d,0,
4418  Comm_pt->mpi_comm());
4419  }
4420  }
4421  }
4422  }
4423 
4425  {
4426  double t_end = TimingHelpers::timer();
4427  oomph_info << "Total time for Mesh::resize_halo_nodes(): "
4428  << t_end-t_start << std::endl;
4429  }
4430 
4431 }
4432 
4433 
4434 //========================================================================
4435 /// Get all the halo data stored in the mesh and add pointers to
4436 /// the data to the map, indexed by global equation number
4437 //=========================================================================
4438 void Mesh::get_all_halo_data(std::map<unsigned,double*> &map_of_halo_data)
4439 {
4440 
4441  //Loop over the map of Halo_nodes
4442  for(std::map<unsigned,Vector<Node*> >::iterator it = Halo_node_pt.begin();
4443  it!=Halo_node_pt.end();++it)
4444  {
4445  //Find the number of nodes
4446  unsigned n_node = (it->second).size();
4447  //Loop over them all
4448  for(unsigned n=0;n<n_node;n++)
4449  {
4450  //Add the Node's values (including any solid values) to
4451  //the map
4452  (it->second)[n]->add_value_pt_to_map(map_of_halo_data);
4453  }
4454  } //End of loop over halo nodes
4455 
4456  //Now loop over all the halo elements and add their internal data
4457  //Loop over the map of Halo_nodes
4458  for(std::map<unsigned,Vector<GeneralisedElement*> >::
4459  iterator it = Root_halo_element_pt.begin();
4460  it!=Root_halo_element_pt.end();++it)
4461  {
4462  //Find the number of root elements
4463  unsigned n_element = (it->second).size();
4464  for(unsigned e=0;e<n_element;e++)
4465  {
4466  GeneralisedElement* el_pt = (it->second)[e];
4467 
4468  // Is it a refineable element?
4469  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(el_pt);
4470  if (ref_el_pt!=0)
4471  {
4472  // Vector of pointers to leaves in tree emanating from
4473  // current root halo element
4474  Vector<Tree*> leaf_pt;
4475  ref_el_pt->tree_pt()->stick_leaves_into_vector(leaf_pt);
4476 
4477  // Loop over leaves and add their objects (the finite elements)
4478  // to vector
4479  unsigned nleaf=leaf_pt.size();
4480  for (unsigned l=0;l<nleaf;l++)
4481  {
4482  leaf_pt[l]->object_pt()->
4483  add_internal_value_pt_to_map(map_of_halo_data);
4484  }
4485  }
4486  else
4487  {
4488  el_pt->add_internal_value_pt_to_map(map_of_halo_data);
4489  }
4490  }
4491  }
4492 
4493  ///Repeat for the external data
4494  for(std::map<unsigned,Vector<Node*> >::iterator it =
4495  External_halo_node_pt.begin();
4496  it!=External_halo_node_pt.end();++it)
4497  {
4498  //Find the number of nodes
4499  unsigned n_node = (it->second).size();
4500  //Loop over them all
4501  for(unsigned n=0;n<n_node;n++)
4502  {
4503  //Add the Node's values (including any solid values) to
4504  //the map
4505  (it->second)[n]->add_value_pt_to_map(map_of_halo_data);
4506  }
4507  } //End of loop over halo nodes
4508 
4509  //Now loop over all the halo elements and add their internal data
4510  //Loop over the map of Halo_nodes
4511  for(std::map<unsigned,Vector<GeneralisedElement*> >::
4512  iterator it = External_halo_element_pt.begin();
4513  it!=External_halo_element_pt.end();++it)
4514  {
4515  //Find the number of root elements
4516  unsigned n_element = (it->second).size();
4517  for(unsigned e=0;e<n_element;e++)
4518  {
4519  (it->second)[e]->add_internal_value_pt_to_map(map_of_halo_data);
4520  }
4521  }
4522 
4523 }
4524 
4525 
4526 
4527 
4528 //========================================================================
4529 /// Get halo node stats for this distributed mesh:
4530 /// Average/max/min number of halo nodes over all processors.
4531 /// \b Careful: Involves MPI Broadcasts and must therefore
4532 /// be called on all processors!
4533 //========================================================================
4534 void Mesh::get_halo_node_stats(double& av_number,
4535  unsigned& max_number,
4536  unsigned& min_number)
4537 {
4538  // Storage for number of processors and current processor
4539  int n_proc=Comm_pt->nproc();
4540  int my_rank=Comm_pt->my_rank();
4541 
4542  // Create vector to hold number of halo nodes
4543  Vector<int> nhalo_nodes(n_proc);
4544 
4545  // Stick own number of halo nodes into appropriate entry
4546  int nhalo_node_local=nhalo_node();
4547 
4548  // Gather information on root processor: First argument group
4549  // specifies what is to be sent (one int from each procssor, indicating
4550  // the number of dofs on it), the second group indicates where
4551  // the results are to be gathered (in rank order) on root processor.
4552  MPI_Gather(&nhalo_node_local,1,MPI_INT,
4553  &nhalo_nodes[0],1, MPI_INT,
4554  0,Comm_pt->mpi_comm());
4555 
4556  // Initialise stats
4557  av_number=0.0;
4558  int max=-1;
4559  int min=1000000000;
4560 
4561  if (my_rank==0)
4562  {
4563  for (int i=0;i<n_proc;i++)
4564  {
4565  av_number+=double(nhalo_nodes[i]);
4566  if (int(nhalo_nodes[i])>max) max=nhalo_nodes[i];
4567  if (int(nhalo_nodes[i])<min) min=nhalo_nodes[i];
4568 
4569  }
4570  av_number/=double(n_proc);
4571  }
4572 
4573  // Now broadcast the result back out
4574  MPI_Bcast(&max,1,MPI_INT,0,Comm_pt->mpi_comm());
4575  MPI_Bcast(&min,1,MPI_INT,0,Comm_pt->mpi_comm());
4576  MPI_Bcast(&av_number,1,MPI_DOUBLE,0,Comm_pt->mpi_comm());
4577 
4578  max_number=max;
4579  min_number=min;
4580 }
4581 
4582 
4583 //========================================================================
4584 /// Get haloed node stats for this distributed mesh:
4585 /// Average/max/min number of haloed nodes over all processors.
4586 /// \b Careful: Involves MPI Broadcasts and must therefore
4587 /// be called on all processors!
4588 //========================================================================
4589  void Mesh::get_haloed_node_stats(double& av_number,
4590  unsigned& max_number,
4591  unsigned& min_number)
4592 {
4593  // Storage for number of processors and current processor
4594  int n_proc=Comm_pt->nproc();
4595  int my_rank=Comm_pt->my_rank();
4596 
4597  // Create vector to hold number of haloed nodes
4598  Vector<int> nhaloed_nodes(n_proc);
4599 
4600  // Stick own number of haloed nodes into appropriate entry
4601  int nhaloed_node_local=nhaloed_node();
4602 
4603  // Gather information on root processor: First argument group
4604  // specifies what is to be sent (one int from each procssor, indicating
4605  // the number of dofs on it), the second group indicates where
4606  // the results are to be gathered (in rank order) on root processor.
4607  MPI_Gather(&nhaloed_node_local,1,MPI_INT,
4608  &nhaloed_nodes[0],1, MPI_INT,
4609  0,Comm_pt->mpi_comm());
4610 
4611  // Initialise stats
4612  av_number=0.0;
4613  int max=-1;
4614  int min=1000000000;
4615 
4616  if (my_rank==0)
4617  {
4618  for (int i=0;i<n_proc;i++)
4619  {
4620  av_number+=double(nhaloed_nodes[i]);
4621  if (int(nhaloed_nodes[i])>max) max=nhaloed_nodes[i];
4622  if (int(nhaloed_nodes[i])<min) min=nhaloed_nodes[i];
4623 
4624  }
4625  av_number/=double(n_proc);
4626  }
4627 
4628  // Now broadcast the result back out
4629  MPI_Bcast(&max,1,MPI_INT,0,Comm_pt->mpi_comm());
4630  MPI_Bcast(&min,1,MPI_INT,0,Comm_pt->mpi_comm());
4631  MPI_Bcast(&av_number,1,MPI_DOUBLE,0,Comm_pt->mpi_comm());
4632 
4633  max_number=max;
4634  min_number=min;
4635 }
4636 
4637 //========================================================================
4638 /// Distribute the mesh. Add to vector of deleted elements.
4639 //========================================================================
4641  const Vector<unsigned>& element_domain,
4642  Vector<GeneralisedElement*>& deleted_element_pt,
4643  DocInfo& doc_info,
4644  const bool& report_stats,
4645  const bool& overrule_keep_as_halo_element_status)
4646  {
4647 
4648  // Store communicator
4649  Comm_pt=comm_pt;
4650 
4651  // Storage for number of processors and current processor
4652  int n_proc=comm_pt->nproc();
4653  int my_rank=comm_pt->my_rank();
4654 
4655  // Storage for number of elements and number of nodes on this mesh
4656  unsigned nelem=this->nelement();
4657  unsigned nnod=this->nnode();
4658 
4659  std::ostringstream filename;
4660 
4661  // Doc the partitioning (only on processor 0)
4662  //-------------------------------------------
4663  if (doc_info.is_doc_enabled())
4664  {
4665  if (my_rank==0)
4666  {
4667  // Open files for doc of element partitioning
4668  Vector<std::ofstream*> domain_file(n_proc);
4669  for (int d=0;d<n_proc;d++)
4670  {
4671  // Note: doc_info.number() was set in Problem::distribute(...) to
4672  // reflect the submesh number
4673  //Clear the filename
4674  filename.str("");
4675  filename << doc_info.directory() << "/domain"
4676  << d << "-" << doc_info.number() << ".dat";
4677  domain_file[d]=new std::ofstream(filename.str().c_str());
4678  }
4679 
4680  // Doc
4681  for (unsigned e=0;e<nelem;e++)
4682  {
4683  //If we can't cast to a finite element, we can't output because
4684  //there is no output function
4685  FiniteElement* f_el_pt =
4686  dynamic_cast<FiniteElement*>(this->element_pt(e));
4687  if(f_el_pt!=0)
4688  {
4689  f_el_pt->output(*domain_file[element_domain[e]],5);
4690  }
4691  }
4692 
4693  for (int d=0;d<n_proc;d++)
4694  {
4695  domain_file[d]->close();
4696  delete domain_file[d];
4697  domain_file[d]=0;
4698  }
4699  }
4700  }
4701 
4702  // Loop over all elements, associate all
4703  //--------------------------------------
4704  // nodes with the highest-numbered processor and record all
4705  //---------------------------------------------------------
4706  // processors the node is associated with
4707  //---------------------------------------
4708 
4709  // Storage for processors in charge and processors associated with data
4710  std::map<Data*,std::set<unsigned> > processors_associated_with_data;
4711  std::map<Data*,unsigned> processor_in_charge;
4712 
4713  // For all nodes set the processor in charge to zero
4714  for (unsigned j=0;j<nnod;j++)
4715  {
4716  Node* nod_pt=this->node_pt(j);
4717  processor_in_charge[nod_pt]=0;
4718  }
4719 
4720  // Loop over elements
4721  for (unsigned e=0;e<nelem;e++)
4722  {
4723  // Get an element and its domain
4724  FiniteElement* el_pt = dynamic_cast<FiniteElement*>(this->element_pt(e));
4725  //Element will only have nodes if it is a finite element
4726  if(el_pt!=0)
4727  {
4728  unsigned el_domain=element_domain[e];
4729 
4730  // Associate nodes with highest numbered processor
4731  unsigned nnod=el_pt->nnode();
4732  for (unsigned j=0;j<nnod;j++)
4733  {
4734  Node* nod_pt=el_pt->node_pt(j);
4735 
4736  // processor in charge was initialised to 0 above
4737  if (el_domain>processor_in_charge[nod_pt])
4738  {
4739  processor_in_charge[nod_pt]=el_domain;
4740  }
4741  processors_associated_with_data[nod_pt].insert(el_domain);
4742  }
4743  }
4744  }
4745 
4746  // Doc the partitioning (only on processor 0)
4747  //-------------------------------------------
4748  if (doc_info.is_doc_enabled())
4749  {
4750  if (my_rank==0)
4751  {
4752  // Open files for doc of node partitioning
4753  Vector<std::ofstream*> node_file(n_proc);
4754  for (int d=0;d<n_proc;d++)
4755  {
4756  // Note: doc_info.number() was set in Problem::distribute(...) to
4757  // reflect the submesh number...
4758  //Clear the filename
4759  filename.str("");
4760  filename << doc_info.directory() << "/node"
4761  << d << "-" << doc_info.number() << ".dat";
4762  node_file[d]=new std::ofstream(filename.str().c_str());
4763  }
4764 
4765  // Doc
4766  for (unsigned j=0;j<nnod;j++)
4767  {
4768  Node* nod_pt=this->node_pt(j);
4769  const unsigned n_dim = nod_pt->ndim();
4770  for(unsigned i=0;i<n_dim;i++)
4771  {
4772  *node_file[processor_in_charge[nod_pt]] << nod_pt->x(i) << " ";
4773  }
4774  *node_file[processor_in_charge[nod_pt]] << "\n";
4775  }
4776  for (int d=0;d<n_proc;d++)
4777  {
4778  node_file[d]->close();
4779  delete node_file[d];
4780  node_file[d]=0;
4781  }
4782  }
4783  }
4784 
4785  // Declare all nodes as obsolete. We'll
4786  // change this setting for all nodes that must be retained
4787  // further down
4788  for (unsigned j=0;j<nnod;j++)
4789  {
4790  this->node_pt(j)->set_obsolete();
4791  }
4792 
4793 
4794  // Backup old mesh data and flush mesh
4795  //-------------------------------------
4796 
4797  // Backup pointers to elements in this mesh
4798  Vector<GeneralisedElement*> backed_up_el_pt(nelem);
4799  for (unsigned e=0;e<nelem;e++)
4800  {
4801  backed_up_el_pt[e]=this->element_pt(e);
4802  }
4803 
4804  // Info. used to re-generate the boundary element scheme after the
4805  // deletion of elements not belonging to the current processor)
4806 
4807  // Get the number of boundary elements before the deletion of non
4808  // retained elements
4809  const unsigned tmp_nboundary = this->nboundary();
4810  Vector<unsigned> ntmp_boundary_elements(tmp_nboundary);
4811  // In case that there are regions, also get the number of boundary
4812  // elements in each region
4813  Vector<Vector<unsigned> > ntmp_boundary_elements_in_region(tmp_nboundary);
4814  // Also get the finite element version of the elements and back them
4815  // up
4816  Vector<FiniteElement*> backed_up_f_el_pt(nelem);
4817 
4818  // Only do this if the mesh is a TriangleMeshBase
4819  TriangleMeshBase* triangle_mesh_pt = dynamic_cast<TriangleMeshBase*>(this);
4820  bool is_a_triangle_mesh_base_mesh = false;
4821  if (triangle_mesh_pt!=0)
4822  {
4823  // Set the flag to indicate we are working with a TriangleMeshBase
4824  // mesh
4825  is_a_triangle_mesh_base_mesh = true;
4826 
4827  // Are there regions?
4828  const unsigned n_regions = triangle_mesh_pt->nregion();
4829 
4830  // Loop over the boundaries
4831  for (unsigned ib = 0; ib < tmp_nboundary; ib++)
4832  {
4833  // Get the number of boundary elements
4834  ntmp_boundary_elements[ib] = this->nboundary_element(ib);
4835 
4836  // Resize the container
4837  ntmp_boundary_elements_in_region[ib].resize(n_regions);
4838 
4839  // Loop over the regions
4840  for (unsigned rr = 0 ; rr < n_regions; rr++)
4841  {
4842  // Get the region id
4843  const unsigned region_id =
4844  static_cast<unsigned>(triangle_mesh_pt->region_attribute(rr));
4845 
4846  // Store the number of element in the region (notice we are
4847  // using the region index not the region id to refer to the
4848  // region)
4849  ntmp_boundary_elements_in_region[ib][rr] =
4850  triangle_mesh_pt->nboundary_element_in_region(ib, region_id);
4851 
4852  } // for (rr < n_regions)
4853 
4854  } // for (ib < tmp_nboundary)
4855 
4856  for (unsigned e=0;e<nelem;e++)
4857  {
4858  // Get the finite element version of the elements and back them
4859  // up
4860  backed_up_f_el_pt[e] = this->finite_element_pt(e);
4861  }
4862 
4863  } // if (triangle_mesh_pt!=0)
4864 
4865  // Flush the mesh storage
4866  this->flush_element_storage();
4867 
4868  // Delete any storage of external elements and nodes
4870 
4871  // Boolean to indicate which element is to be retained
4872  std::vector<bool> element_retained(nelem,false);
4873 
4874  // Storage for element numbers of root halo elements that will be
4875  // retained on current processor: root_halo_element[p][j]
4876  // stores the element number (in the order in which the elements are stored
4877  // in backed_up_el_pt) of the j-th root halo element with processor
4878  // p.
4879  Vector<Vector<int> > root_halo_element(n_proc);
4880 
4881  // Dummy entry to make sure we always have something to send
4882  for (int p=0;p<n_proc;p++)
4883  {
4884  root_halo_element[p].push_back(-1);
4885  }
4886 
4887  // Determine which elements are going to end up on which processor
4888  //----------------------------------------------------------------
4889  unsigned number_of_retained_elements=0;
4890 
4891  // Loop over all backed up elements
4892  nelem=backed_up_el_pt.size();
4893  for (unsigned e=0;e<nelem;e++)
4894  {
4895  // Get element and its domain
4896  GeneralisedElement* el_pt=backed_up_el_pt[e];
4897  unsigned el_domain=element_domain[e];
4898 
4899  // If element is located on current processor add it back to the mesh
4900  if (el_domain==unsigned(my_rank))
4901  {
4902  // Add element to current processor
4903  element_retained[e]=true;
4904  number_of_retained_elements++;
4905  }
4906  // Otherwise we may still need it if it's a halo element:
4907  else
4908  {
4909  // If this current mesh has been told to keep all elements as halos,
4910  // OR the element itself knows that it must be kept then
4911  // keep it
4912  if ((this->Keep_all_elements_as_halos) ||
4913  (el_pt->must_be_kept_as_halo()))
4914  {
4915  if (!overrule_keep_as_halo_element_status)
4916  {
4917  // Add as root halo element whose non-halo counterpart is
4918  // located on processor el_domain
4919  if (!element_retained[e])
4920  {
4921  root_halo_element[el_domain].push_back(e);
4922  element_retained[e]=true;
4923  number_of_retained_elements++;
4924  }
4925  }
4926  }
4927  //Otherwise: Is one of the nodes associated with the current processor?
4928  else
4929  {
4930  //Can only have nodes if this is a finite element
4931  FiniteElement* finite_el_pt = dynamic_cast<FiniteElement*>(el_pt);
4932  if(finite_el_pt!=0)
4933  {
4934  unsigned n_node = finite_el_pt->nnode();
4935  for (unsigned n=0;n<n_node;n++)
4936  {
4937  Node* nod_pt=finite_el_pt->node_pt(n);
4938 
4939  // Keep element? (use stl find?)
4940  unsigned keep_it=false;
4941  for (std::set<unsigned>::iterator
4942  it=processors_associated_with_data[nod_pt].begin();
4943  it!=processors_associated_with_data[nod_pt].end();
4944  it++)
4945  {
4946  if (*it==unsigned(my_rank))
4947  {
4948  keep_it=true;
4949  //Break out of the loop over processors
4950  break;
4951  }
4952  }
4953 
4954  // Add a root halo element either if keep_it=true
4955  if (keep_it)
4956  {
4957  // Add as root halo element whose non-halo counterpart is
4958  // located on processor el_domain
4959  if (!element_retained[e])
4960  {
4961  root_halo_element[el_domain].push_back(e);
4962  element_retained[e]=true;
4963  number_of_retained_elements++;
4964  }
4965  //Now break out of loop over nodes
4966  break;
4967  }
4968  }
4969  }
4970  } //End of testing for halo by virtue of shared nodes
4971  }//End of halo element conditions
4972  } //end of loop over elements
4973 
4974  // First check that the number of elements is greater than zero, when
4975  // working with submeshes it may be possible that some of them have
4976  // no elements (face element meshes) since those may have been
4977  // deleted in "Problem::actions_before_distribute()"
4978  if (nelem > 0)
4979  {
4980  // Check that we are working with a TriangleMeshBase mesh, if
4981  // that is the case then we need to create shared boundaries
4982  if (is_a_triangle_mesh_base_mesh)
4983  {
4984  // Creation of shared boundaries
4985  // ------------------------------
4986  // All processors share the same boundary id for the created
4987  // shared boundary. We need all the elements on all processors,
4988  // that is why this step is performed before the deletion of the
4989  // elements not associated to the current processor.
4990  // N.B.: This applies only to unstructured meshes
4991  this->create_shared_boundaries(comm_pt, element_domain,
4992  backed_up_el_pt,
4993  backed_up_f_el_pt,
4994  processors_associated_with_data,
4995  overrule_keep_as_halo_element_status);
4996  } // if (is_a_triangle_mesh_base_mesh)
4997  } // if (nelem > 0)
4998 
4999  // NOTE: No need to add additional layer of halo elements.
5000  // Procedure for removing "overlooked" halo nodes in
5001  // deals classify_halo_and_haloed_nodes() deals
5002  // with the problem addressed here. [code that dealt with this
5003  // problem at distribution stage has been removed]
5004 
5005  // Store the finite element pointer version of the elements that are
5006  // about to be deleted, used to reset the boundary elements info
5007  Vector<FiniteElement*> deleted_f_element_pt;
5008 
5009  // Copy the elements associated with the actual
5010  // current processor into its own permanent storage.
5011  // Do it in the order in which the elements appeared originally
5012  nelem=backed_up_el_pt.size();
5013  for (unsigned e=0;e<nelem;e++)
5014  {
5015  GeneralisedElement* el_pt=backed_up_el_pt[e];
5016  if (element_retained[e])
5017  {
5018  this->add_element_pt(el_pt);
5019  }
5020  else
5021  {
5022  // Flush the object attached to the tree for this element?
5023  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(el_pt);
5024  if (ref_el_pt!=0)
5025  {
5026  ref_el_pt->tree_pt()->flush_object();
5027  }
5028 
5029 
5030  // Store pointer to the element that's about to be deleted.
5031 
5032  // Only for structured meshes since this "deleted_element_pt"
5033  // vector is used in the "problem" class to set null pointer to
5034  // the deleted elements in the Base_mesh_element_pt structure
5035  if (!is_a_triangle_mesh_base_mesh)
5036  {
5037  deleted_element_pt.push_back(el_pt);
5038  } // if (!is_a_triangle_mesh_base_mesh)
5039 
5040  if (is_a_triangle_mesh_base_mesh)
5041  {
5042  // Store pointer to the finite element that's about to be deleted
5043  deleted_f_element_pt.push_back(backed_up_f_el_pt[e]);
5044  }
5045 
5046  // Delete the element
5047  delete el_pt;
5048  }
5049  }
5050 
5051  // Copy the root halo elements associated with the actual
5052  // current processor into its own permanent storage; the order
5053  // here is somewhat random but we compensate for that by
5054  // ensuring that the corresponding haloed elements are
5055  // added in the same order below
5056 #ifdef PARANOID
5057  std::map<unsigned,bool> done;
5058 #endif
5059  for (int d=0;d<n_proc;d++)
5060  {
5061  nelem=root_halo_element[d].size();
5062  for (unsigned e=0;e<nelem;e++)
5063  {
5064  int number=root_halo_element[d][e];
5065  if (number>=0)
5066  {
5067 #ifdef PARANOID
5068  if (done[number])
5069  {
5070  std::ostringstream error_message;
5071  error_message
5072  << "Have already added element " << number
5073  << " as root halo element\n"
5074  << std::endl;
5075  throw OomphLibError(error_message.str(),
5076  OOMPH_CURRENT_FUNCTION,
5077  OOMPH_EXCEPTION_LOCATION);
5078  }
5079  done[number]=true;
5080 #endif
5081  this->add_root_halo_element_pt(d,backed_up_el_pt[number]);
5082  }
5083  }
5084  }
5085 
5086 
5087  // Now get root haloed elements: root_haloed_element[p][j] stores
5088  // the element number (in the order in which the elements are stored
5089  // in backedup_el_pt) of the j-th rooted halo element with processor
5090  // p. On proc my_proc this the same as root_haloed_element[my_proc][j]
5091  // on processor p, so get the information by gatherv operations.
5092  Vector<Vector<unsigned> > root_haloed_element(n_proc);
5093 
5094  // Find out number of root halo elements with each other proc
5095  Vector<int> nhalo(n_proc,0);
5096  Vector<int> nhaloed(n_proc,0);
5097  for (int p=0;p<n_proc;p++)
5098  {
5099  nhalo[p]=root_halo_element[p].size();
5100  }
5101 
5102  // Each processor sends number of halo elements it has with processor
5103  // p to processor p where this information is stored in nhaloed[...]
5104  for (int p=0;p<n_proc;p++)
5105  {
5106  // Gather the p-th entries in nhalo from every processor on
5107  // processor p and store them in nhaloed consecutively
5108  // starting at beginning
5109  MPI_Gather(&nhalo[p], 1, MPI_INT,
5110  &nhaloed[0], 1, MPI_INT,
5111  p,comm_pt->mpi_comm());
5112  }
5113 
5114  // In the big sequence of concatenated root halo elements (enumerated
5115  // individually on the various processors) where do the root halo
5116  // elements from a given processor start? Also figure out how many
5117  // root haloed elements there are in total by summing up their numbers
5118  Vector<int> start_index(n_proc,0);
5119  unsigned total_number_of_root_haloed_elements=0;
5120  for (int i_proc=0; i_proc<n_proc; i_proc++)
5121  {
5122  total_number_of_root_haloed_elements+=nhaloed[i_proc];
5123  if (i_proc!=0)
5124  {
5125  start_index[i_proc]=total_number_of_root_haloed_elements-
5126  nhaloed[i_proc];
5127  }
5128  else
5129  {
5130  start_index[0]=0;
5131  }
5132  }
5133 
5134  // Storage for all root haloed elements from the various processors, one
5135  // after the other, with some padding from negative entries to avoid
5136  // zero length vectors
5137  Vector<int> all_root_haloed_element(total_number_of_root_haloed_elements,
5138  0);
5139 
5140  // Now send the ids of the relevant elements via gatherv
5141  for (int p=0;p<n_proc;p++)
5142  {
5143  // Gather the p-th entries in nhalo from every processor on
5144  // processor p and store them in nhaloed consecutively
5145  // starting at beginning
5146  MPI_Gatherv(&root_halo_element[p][0], // pointer to first entry in vector
5147  // to be gathered on processor p
5148  nhalo[p], // Number of entries to be sent
5149  MPI_INT,
5150  &all_root_haloed_element[0], // Target -- this will store
5151  // the element numbers of
5152  // all root haloed elements
5153  // received from other processors
5154  // in order
5155  &nhaloed[0], // Pointer to the vector containing the lengths
5156  // of the vectors received from elsewhere
5157  &start_index[0], // "offset" for storage of vector received
5158  // from various processors in the global
5159  // concatenated vector
5160  MPI_INT,
5161  p, // processor that gathers the information
5162  comm_pt->mpi_comm());
5163  }
5164 
5165 
5166  // Determine root haloed elements
5167  //-------------------------------
5168 
5169  // Loop over all other processors
5170  unsigned count=0;
5171  for (int d=0;d<n_proc;d++)
5172  {
5173 
5174 #ifdef PARANOID
5175  std::map<unsigned,bool> done;
5176 #endif
5177 
5178  // Loop over root haloed elements with specified processor
5179  unsigned n=nhaloed[d];
5180  for (unsigned e=0;e<n;e++)
5181  {
5182 
5183  int number=all_root_haloed_element[count];
5184  count++;
5185 
5186  // Ignore padded -1s that were only inserted to avoid
5187  // zero sized vectors
5188  if (number>=0)
5189  {
5190  // Get pointer to element
5191  GeneralisedElement* el_pt=backed_up_el_pt[number];
5192 
5193  // Halo elements can't be haloed themselves
5194  if (!el_pt->is_halo())
5195  {
5196 
5197 #ifdef PARANOID
5198  if (done[number])
5199  {
5200  std::ostringstream error_message;
5201  error_message
5202  << "Have already added element " << number
5203  << " as root haloed element\n"
5204  << std::endl;
5205  throw OomphLibError(error_message.str(),
5206  OOMPH_CURRENT_FUNCTION,
5207  OOMPH_EXCEPTION_LOCATION);
5208  }
5209  done[number]=true;
5210 #endif
5211 
5212  // Current element is haloed by other processor
5213  this->add_root_haloed_element_pt(d,el_pt);
5214  }
5215  }
5216  }
5217  }
5218 
5219 
5220  // Doc stats
5221  if (report_stats)
5222  {
5223  oomph_info << "Processor " << my_rank
5224  << " holds " << this->nelement()
5225  << " elements of which " << this->nroot_halo_element()
5226  << " are root halo elements \n while "
5227  << this->nroot_haloed_element()
5228  << " are root haloed elements" << std::endl;
5229  }
5230 
5231 
5232  // Loop over all retained elements and mark their nodes
5233  //-----------------------------------------------------
5234  // as to be retained too (some double counting going on here)
5235  //-----------------------------------------------------------
5236  nelem=this->nelement();
5237  for (unsigned e=0;e<nelem;e++)
5238  {
5239  FiniteElement* f_el_pt= dynamic_cast<FiniteElement*>(this->element_pt(e));
5240 
5241  //If we have a finite element
5242  if(f_el_pt!=0)
5243  {
5244  // Loop over nodes
5245  unsigned nnod=f_el_pt->nnode();
5246  for (unsigned j=0;j<nnod;j++)
5247  {
5248  Node* nod_pt=f_el_pt->node_pt(j);
5249  nod_pt->set_non_obsolete();
5250  }
5251  }
5252  }
5253 
5254 
5255  // Now remove the pruned nodes
5256  this->prune_dead_nodes();
5257 
5258 
5259 #ifdef OOMPH_HAS_TRIANGLE_LIB
5260  if (is_a_triangle_mesh_base_mesh)
5261  {
5262  triangle_mesh_pt->
5263  reset_boundary_element_info(ntmp_boundary_elements,
5264  ntmp_boundary_elements_in_region,
5265  deleted_f_element_pt);
5266  } // if (tri_mesh_pt!=0)
5267  else
5268  {
5269 #endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
5270 
5271  // Temporarly set the mesh as distributed
5273 
5274 #ifdef OOMPH_HAS_TRIANGLE_LIB
5275  }
5276 #endif
5277 
5278  // Re-setup tree forest. (Call this every time even if
5279  // a (distributed) mesh has no elements on this processor.
5280  // We still need to participate in communication.)
5281  TreeBasedRefineableMeshBase* ref_mesh_pt=
5282  dynamic_cast<TreeBasedRefineableMeshBase*>(this);
5283  if (ref_mesh_pt!=0)
5284  {
5285  ref_mesh_pt->setup_tree_forest();
5286  }
5287 
5288  // Classify nodes
5289  classify_halo_and_haloed_nodes(doc_info,report_stats);
5290 
5291  // Doc?
5292  //-----
5293  if (doc_info.is_doc_enabled())
5294  {
5295  doc_mesh_distribution(doc_info);
5296  }
5297 
5298 }
5299 
5300 
5301 //========================================================================
5302 /// (Irreversibly) prune halo(ed) elements and nodes, usually
5303 /// after another round of refinement, to get rid of
5304 /// excessively wide halo layers. Note that the current
5305 /// mesh will be now regarded as the base mesh and no unrefinement
5306 /// relative to it will be possible once this function
5307 /// has been called.
5308 //========================================================================
5310  deleted_element_pt,
5311  DocInfo& doc_info,
5312  const bool& report_stats)
5313 {
5314 
5315  //MemoryUsage::doc_memory_usage("at start of mesh-level prunes");
5316 
5317 
5318 #ifdef OOMPH_HAS_MPI
5319  // Delete any external element storage before performing the redistribution
5320  // (in particular, external halo nodes that are on mesh boundaries)
5322 #endif
5323 
5324  TreeBasedRefineableMeshBase* ref_mesh_pt=
5325  dynamic_cast<TreeBasedRefineableMeshBase*>(this);
5326  if (ref_mesh_pt!=0)
5327  {
5328 
5329  // Storage for number of processors and current processor
5330  int n_proc=Comm_pt->nproc();
5331  int my_rank=Comm_pt->my_rank();
5332 
5333  // Doc stats
5334  if (report_stats)
5335  {
5336  oomph_info << "\n----------------------------------------------------\n";
5337  oomph_info << "Before pruning: Processor " << my_rank
5338  << " holds " << this->nelement()
5339  << " elements of which " << this->nroot_halo_element()
5340  << " are root halo elements \n while "
5341  << this->nroot_haloed_element()
5342  << " are root haloed elements" << std::endl;
5343 
5344  // Report total number of halo(ed) and shared nodes for this process
5345  oomph_info << "Before pruning: Processor " << my_rank
5346  << " holds " << this->nnode()
5347  << " nodes of which " << this->nhalo_node()
5348  << " are halo nodes \n while " << this->nhaloed_node()
5349  << " are haloed nodes, and " << this->nshared_node()
5350  << " are shared nodes." << std::endl;
5351 
5352  // Report number of halo(ed) and shared nodes with each domain
5353  // from the current process
5354  for (int iproc=0;iproc<n_proc;iproc++)
5355  {
5356  oomph_info << "Before pruning: With process " << iproc << ", there are "
5357  << this->nhalo_node(iproc) << " halo nodes, and " << std::endl
5358  << this->nhaloed_node(iproc) << " haloed nodes, and "
5359  << this->nshared_node(iproc) << " shared nodes" << std::endl;
5360  }
5361  oomph_info << "----------------------------------------------------\n\n";
5362  }
5363 
5364 
5365  double t_start = 0.0;
5367  {
5368  t_start=TimingHelpers::timer();
5369  }
5370 
5371  // Declare all nodes as obsolete. We'll
5372  // change this setting for all nodes that must be retained
5373  // further down
5374  unsigned nnod=this->nnode();
5375  for (unsigned j=0;j<nnod;j++)
5376  {
5377  this->node_pt(j)->set_obsolete();
5378  }
5379 
5380  // Backup pointers to elements in this mesh (they must be finite elements
5381  // beacuse it's a refineable mesh)
5382  unsigned nelem=this->nelement();
5383  Vector<FiniteElement*> backed_up_el_pt(nelem);
5384  std::map<RefineableElement*,bool> keep_element;
5385  for (unsigned e=0;e<nelem;e++)
5386  {
5387  backed_up_el_pt[e]=this->finite_element_pt(e);
5388  }
5389 
5390  //MemoryUsage::doc_memory_usage("after backed up elements");
5391 
5392  // Get the min and max refinement level, and current refinement pattern
5393  unsigned min_ref=0;
5394  unsigned max_ref=0;
5395 
5396  // Get min and max refinement level
5397  unsigned local_min_ref=0;
5398  unsigned local_max_ref=0;
5399  ref_mesh_pt->get_refinement_levels(local_min_ref,local_max_ref);
5400 
5401  // Reconcile between processors: If (e.g. following distribution/pruning)
5402  // the mesh has no elements on this processor) then ignore its
5403  // contribution to the poll of max/min refinement levels
5404  int int_local_min_ref=local_min_ref;
5405  int int_local_max_ref=local_max_ref;
5406 
5407  if (nelem==0)
5408  {
5409  int_local_min_ref=INT_MAX;
5410  int_local_max_ref=INT_MIN;
5411  }
5412 
5413  int int_min_ref=0;
5414  int int_max_ref=0;
5415 
5416  MPI_Allreduce(&int_local_min_ref,&int_min_ref,1,
5417  MPI_INT,MPI_MIN,
5418  Comm_pt->mpi_comm());
5419  MPI_Allreduce(&int_local_max_ref,&int_max_ref,1,
5420  MPI_INT,MPI_MAX,
5421  Comm_pt->mpi_comm());
5422 
5423  min_ref=unsigned(int_min_ref);
5424  max_ref=unsigned(int_max_ref);
5425 
5426  // Get refinement pattern
5427  Vector<Vector<unsigned> > current_refined;
5428  ref_mesh_pt->get_refinement_pattern(current_refined);
5429 
5430  // get_refinement_pattern refers to the elements at each level
5431  // that were refined when proceeding to the next level
5432  int local_n_ref=current_refined.size();
5433 
5434  // Bypass if no elements
5435  if (nelem==0)
5436  {
5437  local_n_ref=INT_MIN;
5438  }
5439 
5440  // Reconcile between processors
5441  int int_n_ref=0;
5442  MPI_Allreduce(&local_n_ref,&int_n_ref,1,
5443  MPI_INT,MPI_MAX,
5444  Comm_pt->mpi_comm());
5445  unsigned n_ref=int(int_n_ref);
5446 
5447 
5448  double t_end = 0.0;
5450  {
5451  t_end=TimingHelpers::timer();
5452  oomph_info
5453  << "Time for establishing refinement levels in "
5454  << " Mesh::prune_halo_elements_and_nodes() [includes comms]: "
5455  << t_end-t_start << std::endl;
5456  t_start = TimingHelpers::timer();
5457  }
5458 
5459  //MemoryUsage::doc_memory_usage("after establishing refinement levels");
5460 
5461  // Bypass everything until next comms if no elements
5462  if (nelem>0)
5463  {
5464  // Loop over all elements; keep those on the min refinement level
5465  // Need to go back to the level indicated by min_ref
5466  unsigned base_level=n_ref-(max_ref-min_ref);
5467 
5468  // This is the new minimum uniform refinement level
5469  // relative to total unrefined base mesh
5470  ref_mesh_pt->uniform_refinement_level_when_pruned()=min_ref;
5471 
5472  // Get the elements at the specified "base" refinement level
5473  Vector<RefineableElement*> base_level_elements_pt;
5474  ref_mesh_pt->get_elements_at_refinement_level(base_level,
5475  base_level_elements_pt);
5476 
5477  // Loop over the elements at this level
5478  unsigned n_base_el=base_level_elements_pt.size();
5479  for (unsigned e=0;e<n_base_el;e++)
5480  {
5481  // Extract correct element...
5482  RefineableElement* ref_el_pt=base_level_elements_pt[e];
5483 
5484 
5485  // Check it exists
5486  if (ref_el_pt!=0)
5487  {
5488  // Keep all non-halo elements, remove excess halos
5489  if ((!ref_el_pt->is_halo())||(ref_el_pt->must_be_kept_as_halo()))
5490  {
5491  keep_element[ref_el_pt]=true;
5492 
5493  // Loop over this non-halo element's nodes and retain them too
5494  unsigned nnod=ref_el_pt->nnode();
5495  for (unsigned j=0;j<nnod;j++)
5496  {
5497  ref_el_pt->node_pt(j)->set_non_obsolete();
5498  }
5499  }
5500  }
5501 
5502  } // end loop over base level elements
5503  }
5504 
5505  // Synchronise refinement level when pruned over all meshes even if they
5506  // were empty (in which case the uniform refinement level is still zero
5507  // so go for max
5508  unsigned n_unif=0;
5509  unsigned n_unif_local=ref_mesh_pt->uniform_refinement_level_when_pruned();
5510  MPI_Allreduce(&n_unif_local,&n_unif,1,
5511  MPI_UNSIGNED,MPI_MAX,
5512  Comm_pt->mpi_comm());
5513  ref_mesh_pt->uniform_refinement_level_when_pruned()=n_unif;
5514 
5515 
5516  t_end = 0.0;
5518  {
5519  t_end=TimingHelpers::timer();
5520  oomph_info
5521  << "Time for synchronising refinement levels in "
5522  << " Mesh::prune_halo_elements_and_nodes() [includes comms]: "
5523  << t_end-t_start << std::endl;
5524  t_start = TimingHelpers::timer();
5525  }
5526 
5527 
5528 
5529  //MemoryUsage::doc_memory_usage("after synchronising refinement levels");
5530 
5531 
5532  // Now work on which "root" halo elements to keep at this level
5533  // Can't use the current set directly; however,
5534  // we know the refinement level of the current halo, so
5535  // it is possible to go from that backwards to find the "father
5536  // halo element" necessary to complete this step
5537 
5538  // Temp map of vectors holding the pointers to the root halo elements
5539  std::map<unsigned, Vector<FiniteElement*> > tmp_root_halo_element_pt;
5540 
5541  // Temp map of vectors holding the pointers to the root haloed elements
5542  std::map<unsigned, Vector<FiniteElement*> > tmp_root_haloed_element_pt;
5543 
5544  // Make sure we only push back each element once
5545  std::map<unsigned,std::map<FiniteElement*,bool> >
5546  tmp_root_halo_element_is_retained;
5547 
5548  // Map to store if a halo element survives
5549  std::map<FiniteElement*,bool> halo_element_is_retained;
5550 
5551  for (int domain=0;domain<n_proc;domain++)
5552  {
5553  // Get vector of halo elements with processor domain by copy operation
5554  Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(domain));
5555 
5556  // Loop over halo elements associated with this adjacent domain
5557  unsigned nelem=halo_elem_pt.size();
5558  for (unsigned e=0;e<nelem;e++)
5559  {
5560  // Get element
5561  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>
5562  (halo_elem_pt[e]);
5563 
5564  // An element should only be kept if its refinement
5565  // level is the same as the minimum refinement level
5566  unsigned halo_el_level=ref_el_pt->refinement_level();
5567 
5568  RefineableElement* el_pt=0;
5569  if (halo_el_level==min_ref)
5570  {
5571  // Already at the correct level
5572  el_pt=ref_el_pt;
5573  }
5574  else
5575  {
5576  // Need to go up the tree to the father element at min_ref
5577  RefineableElement* father_el_pt;
5578  ref_el_pt->get_father_at_refinement_level(min_ref,father_el_pt);
5579  el_pt=father_el_pt;
5580  }
5581 
5582  // Now loop over nodes in the halo element and retain it as
5583  // halo element if any of it's nodes are shared with one of the
5584  // non-halo-elements that we've already identified earlier -- these
5585  // were set to non-obsolete above.
5586  unsigned nnod=el_pt->nnode();
5587  for (unsigned j=0;j<nnod;j++)
5588  {
5589  // Node has been reclaimed by one of the non-halo elements above
5590  Node* nod_pt=el_pt->node_pt(j);
5591  if (!nod_pt->is_obsolete())
5592  {
5593  // Keep element and add it to preliminary storage for
5594  // halo elements associated with current neighbouring domain
5595  if (!tmp_root_halo_element_is_retained[domain][el_pt])
5596  {
5597  tmp_root_halo_element_pt[domain].push_back(el_pt);
5598  tmp_root_halo_element_is_retained[domain][el_pt]=true;
5599  }
5600  keep_element[el_pt]=true;
5601  halo_element_is_retained[el_pt]=true;
5602  break;
5603  }
5604  }
5605  }
5606  }
5607 
5609  {
5610  t_end=TimingHelpers::timer();
5611  oomph_info
5612  << "Time for setup of retention pattern in "
5613  << " Mesh::prune_halo_elements_and_nodes(): "
5614  << t_end-t_start << std::endl;
5615  t_start = TimingHelpers::timer();
5616  }
5617 
5618  //MemoryUsage::doc_memory_usage("after setup of retention pattern");
5619 
5620  // Make sure everybody finishes this part
5621  MPI_Barrier(Comm_pt->mpi_comm());
5622 
5623  // Now all processors have decided (independently) which of their
5624  // (to-be root) halo elements they wish to retain. Now we need to figure out
5625  // which of their elements are haloed and add them in the appropriate
5626  // order into the haloed element scheme. For this we exploit that
5627  // the halo and haloed elements are accessed in the same order on
5628  // all processors!
5629 
5630  // Identify haloed elements on domain d
5631  for (int d=0;d<n_proc;d++)
5632  {
5633  // Loop over domains that halo this domain
5634  for (int dd=0;dd<n_proc;dd++)
5635  {
5636  // Dont't talk to yourself
5637  if (d!=dd)
5638  {
5639  // If we're identifying my haloed elements:
5640  if (d==my_rank)
5641  {
5642  // Get vector all elements that are currently haloed by domain dd
5644  haloed_elem_pt(this->haloed_element_pt(dd));
5645 
5646  // Create a vector of ints to indicate if the halo element
5647  // on processor dd processor was kept
5648  unsigned nelem=haloed_elem_pt.size();
5649  Vector<int> halo_kept(nelem);
5650 
5651  // Receive this vector from processor dd
5652  if (nelem!=0)
5653  {
5654  MPI_Status status;
5655  MPI_Recv(&halo_kept[0],nelem,MPI_INT,dd,0,Comm_pt->mpi_comm(),
5656  &status);
5657 
5658  // Classify haloed element accordingly
5659  for (unsigned e=0;e<nelem;e++)
5660  {
5661  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>
5662  (haloed_elem_pt[e]);
5663 
5664  // An element should only be kept if its refinement
5665  // level is the same as the minimum refinement level
5666  unsigned haloed_el_level=ref_el_pt->refinement_level();
5667 
5668  // Go up the tree to the correct level
5669  RefineableElement* el_pt;
5670 
5671  if (haloed_el_level==min_ref)
5672  {
5673  // Already at the correct level
5674  el_pt=ref_el_pt;
5675  }
5676  else
5677  {
5678  // Need to go up the tree to the father element at min_ref
5679  RefineableElement* father_el_pt;
5681  (min_ref,father_el_pt);
5682  el_pt=father_el_pt;
5683  }
5684 
5685  if (halo_kept[e]==1)
5686  {
5687  // I am being haloed by processor dd
5688  // Only keep it if it's not already in the storage
5689  bool already_root_haloed=false;
5690  unsigned n_root_haloed=tmp_root_haloed_element_pt[dd].size();
5691  for (unsigned e_root=0;e_root<n_root_haloed;e_root++)
5692  {
5693  if (el_pt==tmp_root_haloed_element_pt[dd][e_root])
5694  {
5695  already_root_haloed=true;
5696  break;
5697  }
5698  }
5699  if (!already_root_haloed)
5700  {
5701  tmp_root_haloed_element_pt[dd].push_back(el_pt);
5702  }
5703  }
5704  }
5705  }
5706  }
5707  else
5708  {
5709  // If we're dealing with my halo elements:
5710  if (dd==my_rank)
5711  {
5712  // Find (current) halo elements on processor dd whose non-halo is
5713  // on processor d
5715  halo_elem_pt(this->halo_element_pt(d));
5716 
5717  // Create a vector of ints to indicate if the halo
5718  // element was kept
5719  unsigned nelem=halo_elem_pt.size();
5720  Vector<int> halo_kept(nelem,0);
5721  for (unsigned e=0;e<nelem;e++)
5722  {
5723  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>
5724  (halo_elem_pt[e]);
5725 
5726  // An element should only be kept if its refinement
5727  // level is the same as the minimum refinement level
5728  unsigned halo_el_level=ref_el_pt->refinement_level();
5729 
5730  // Go up the tree to the correct level
5731  RefineableElement* el_pt;
5732  if (halo_el_level==min_ref)
5733  {
5734  // Already at the correct level
5735  el_pt=ref_el_pt;
5736  }
5737  else
5738  {
5739  // Need to go up the tree to the father element at min_ref
5740  RefineableElement* father_el_pt;
5742  (min_ref,father_el_pt);
5743  el_pt=father_el_pt;
5744  }
5745 
5746  if (halo_element_is_retained[el_pt])
5747  {
5748  halo_kept[e]=1;
5749  }
5750  }
5751 
5752  // Now send this vector to processor d to tell it which of
5753  // the haloed elements (which are listed in the same order)
5754  // are to be retained as haloed elements.
5755  if (nelem!=0)
5756  {
5757  MPI_Send(&halo_kept[0],nelem,MPI_INT,d,0,Comm_pt->mpi_comm());
5758  }
5759  }
5760  }
5761  }
5762  }
5763  }
5764 
5766  {
5767  t_end=TimingHelpers::timer();
5768  oomph_info
5769  << "Time for pt2pt comms of retention pattern in "
5770  << " Mesh::prune_halo_elements_and_nodes(): "
5771  << t_end-t_start << std::endl;
5772  t_start = TimingHelpers::timer();
5773  }
5774 
5775  //MemoryUsage::doc_memory_usage("after pt2pt comms of retention pattern");
5776 
5777 
5778  // Backup pointers to nodes in this mesh
5779  nnod=this->nnode();
5780  Vector<Node*> backed_up_nod_pt(nnod);
5781  for (unsigned j=0;j<nnod;j++)
5782  {
5783  backed_up_nod_pt[j]=this->node_pt(j);
5784  }
5785 
5786  // Flush the mesh storage
5788 
5789  // Storage for non-active trees/elements that are to be deleted
5790  std::set<Tree*> trees_to_be_deleted_pt;
5791 
5792  // Loop over all backed-up elements
5793  nelem=backed_up_el_pt.size();
5794  for (unsigned e=0;e<nelem;e++)
5795  {
5796  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>
5797  (backed_up_el_pt[e]);
5798 
5799  // Get refinement level
5800  unsigned level=ref_el_pt->refinement_level();
5801 
5802  // Go up the tree to the correct level
5803  RefineableElement* el_pt;
5804 
5805  if (level==min_ref)
5806  {
5807  // Already at the correct level
5808  el_pt=ref_el_pt;
5809  }
5810  else
5811  {
5812  // Need to go up the tree to the father element at min_ref
5813  RefineableElement* father_el_pt;
5815  (min_ref,father_el_pt);
5816  el_pt=father_el_pt;
5817  }
5818 
5819  // If the base element is going to be kept, then add the current element
5820  // to the "new" mesh
5821  if (keep_element[el_pt])
5822  {
5823  this->add_element_pt(ref_el_pt);
5824  }
5825  else
5826  {
5827  // Flush the object attached to the tree for this element?
5828  RefineableElement* my_el_pt=dynamic_cast<RefineableElement*>(ref_el_pt);
5829  if (my_el_pt!=0)
5830  {
5831  my_el_pt->tree_pt()->flush_object();
5832  }
5833 
5834  // Get associated tree root
5835  Tree* tmp_tree_root_pt=my_el_pt->tree_pt()->root_pt();
5836 
5837  // Get all the tree nodes
5838  Vector<Tree*> tmp_all_tree_nodes_pt;
5839  tmp_tree_root_pt->stick_all_tree_nodes_into_vector(
5840  tmp_all_tree_nodes_pt);
5841 
5842  // Loop over all of them and delete associated object/
5843  // and flush any record of it in tree
5844  unsigned n_tree=tmp_all_tree_nodes_pt.size();
5845  for (unsigned j=0;j<n_tree;j++)
5846  {
5847  if (tmp_all_tree_nodes_pt[j]->object_pt()!=0)
5848  {
5849  unsigned lev=tmp_all_tree_nodes_pt[j]->object_pt()->refinement_level();
5850  if (lev<=min_ref)
5851  {
5852  if (!keep_element[tmp_all_tree_nodes_pt[j]->object_pt()])
5853  {
5854  trees_to_be_deleted_pt.insert(tmp_all_tree_nodes_pt[j]);
5855  }
5856  }
5857  }
5858  }
5859 
5860  // Delete the element
5861  deleted_element_pt.push_back(ref_el_pt);
5862  delete ref_el_pt;
5863  }
5864  }
5865 
5866  //MemoryUsage::doc_memory_usage("before deleting superfluous elements");
5867 
5868  // Delete superfluous elements
5869  for (std::set<Tree*>::iterator it=trees_to_be_deleted_pt.begin();
5870  it!=trees_to_be_deleted_pt.end();it++)
5871  {
5872  Tree* tree_pt=(*it);
5873  if (tree_pt->object_pt()!=0)
5874  {
5875  deleted_element_pt.push_back(tree_pt->object_pt());
5876  delete tree_pt->object_pt();
5877  tree_pt->flush_object();
5878  }
5879  }
5880 
5881  //MemoryUsage::doc_memory_usage("after deleting superfluous elements");
5882 
5883 
5884  // Wipe the storage scheme for (root) halo(ed) elements and then re-assign
5885  Root_haloed_element_pt.clear();
5886  Root_halo_element_pt.clear();
5887  for (int domain=0;domain<n_proc;domain++)
5888  {
5889  unsigned nelem=tmp_root_halo_element_pt[domain].size();
5890  for (unsigned e=0;e<nelem;e++)
5891  {
5892  Root_halo_element_pt[domain].push_back(
5893  tmp_root_halo_element_pt[domain][e]);
5894  }
5895 
5896  nelem=tmp_root_haloed_element_pt[domain].size();
5897  for (unsigned e=0;e<nelem;e++)
5898  {
5899  Root_haloed_element_pt[domain].push_back(
5900  tmp_root_haloed_element_pt[domain][e]);
5901  }
5902  }
5903 
5904 // MemoryUsage::doc_memory_usage(
5905 // "after wiping storage scheme for root halo/ed els");
5906 
5907  // Loop over all retained elements at this level and mark their nodes
5908  //-------------------------------------------------------------------
5909  // as to be retained too (some double counting going on here)
5910  //-----------------------------------------------------------
5911  nelem=this->nelement();
5912  for (unsigned e=0;e<nelem;e++)
5913  {
5914  FiniteElement* el_pt=this->finite_element_pt(e);
5915 
5916  // Loop over nodes
5917  unsigned nnod=el_pt->nnode();
5918  for (unsigned j=0;j<nnod;j++)
5919  {
5920  Node* nod_pt=el_pt->node_pt(j);
5921  nod_pt->set_non_obsolete();
5922  }
5923  }
5924 
5925  // Complete rebuild of mesh by adding retained nodes
5926  // Note that they are added in the order in which they
5927  // occured in the original mesh as this guarantees the
5928  // synchronisity between the serialised access to halo
5929  // and haloed nodes from different processors.
5930  nnod=backed_up_nod_pt.size();
5931  for (unsigned j=0;j<nnod;j++)
5932  {
5933  Node* nod_pt=backed_up_nod_pt[j];
5934  if(!nod_pt->is_obsolete())
5935  {
5936  // Not obsolete so add it back to the mesh
5937  this->add_node_pt(nod_pt);
5938  }
5939  }
5940 
5941  // MemoryUsage::doc_memory_usage("after adding nodes back in");
5942 
5943  // Prune and rebuild mesh
5944  //-----------------------
5945 
5946  // Now remove the pruned nodes from the boundary lookup scheme
5947  this->prune_dead_nodes();
5948 
5949  // MemoryUsage::doc_memory_usage("after prune dead nodes");
5950 
5951  // And finally re-setup the boundary lookup scheme for elements
5953 
5954  //MemoryUsage::doc_memory_usage("after setup boundary info");
5955 
5956  // Re-setup tree forest if needed. (Call this every time even if
5957  // a (distributed) mesh has no elements on this processor.
5958  // We still need to participate in communication.)
5959  TreeBasedRefineableMeshBase* ref_mesh_pt=
5960  dynamic_cast<TreeBasedRefineableMeshBase*>(this);
5961  if (ref_mesh_pt!=0)
5962  {
5963  ref_mesh_pt->setup_tree_forest();
5964  }
5965 
5966  //MemoryUsage::doc_memory_usage("after setup tree forest");
5967 
5969  {
5970  t_end=TimingHelpers::timer();
5971  oomph_info
5972  << "Time for local rebuild of mesh from retained elements in "
5973  << " Mesh::prune_halo_elements_and_nodes(): "
5974  << t_end-t_start << std::endl;
5975  t_start = TimingHelpers::timer();
5976  }
5977 
5978  // Classify nodes
5979  classify_halo_and_haloed_nodes(doc_info,report_stats);
5980 
5982  {
5983  t_end=TimingHelpers::timer();
5984  oomph_info
5985  << "Time for Mesh::classify_halo_and_haloed_nodes() "
5986  << " Mesh::prune_halo_elements_and_nodes(): "
5987  << t_end-t_start << std::endl;
5988  t_start = TimingHelpers::timer();
5989  }
5990 
5991  //MemoryUsage::doc_memory_usage("after classify_halo_and_haloed_nodes()");
5992 
5993  // Doc?
5994  //-----
5995  if (doc_info.is_doc_enabled())
5996  {
5997  oomph_info << "Outputting distribution in "
5998  << doc_info.directory() << " "
5999  << doc_info.number() << std::endl;
6000  doc_mesh_distribution(doc_info);
6001  }
6002 
6003  // Finally: Reorder the nodes within the mesh's node vector
6004  // to establish a standard ordering regardless of the sequence
6005  // of mesh refinements -- this is required to allow dump/restart
6006  // on refined meshes
6007  this->reorder_nodes();
6008 
6009 
6011  {
6012  t_end=TimingHelpers::timer();
6013  oomph_info
6014  << "Time for Mesh::reorder_nodes() "
6015  << " Mesh::prune_halo_elements_and_nodes(): "
6016  << t_end-t_start << std::endl;
6017  t_start = TimingHelpers::timer();
6018  }
6019 
6020  //MemoryUsage::doc_memory_usage("after reorder nodes");
6021 
6022  // Doc stats
6023  if (report_stats)
6024  {
6025  oomph_info << "\n----------------------------------------------------\n";
6026  oomph_info << "After pruning: Processor " << my_rank
6027  << " holds " << this->nelement()
6028  << " elements of which " << this->nroot_halo_element()
6029  << " are root halo elements \n while "
6030  << this->nroot_haloed_element()
6031  << " are root haloed elements" << std::endl;
6032 
6033  // Report total number of halo(ed) and shared nodes for this process
6034  oomph_info << "After pruning: Processor " << my_rank
6035  << " holds " << this->nnode()
6036  << " nodes of which " << this->nhalo_node()
6037  << " are halo nodes \n while " << this->nhaloed_node()
6038  << " are haloed nodes, and " << this->nshared_node()
6039  << " are shared nodes." << std::endl;
6040 
6041  // Report number of halo(ed) and shared nodes with each domain
6042  // from the current process
6043  for (int iproc=0;iproc<n_proc;iproc++)
6044  {
6045  oomph_info << "After pruning: With process " << iproc << ", there are "
6046  << this->nhalo_node(iproc) << " halo nodes, and " << std::endl
6047  << this->nhaloed_node(iproc) << " haloed nodes, and "
6048  << this->nshared_node(iproc) << " shared nodes" << std::endl;
6049  }
6050  oomph_info << "----------------------------------------------------\n\n";
6051  }
6052 
6053  }
6054 
6055  //MemoryUsage::doc_memory_usage("end of mesh level prune");
6056 
6057 }
6058 
6059 
6060 
6061 
6062 
6063 
6064 //========================================================================
6065 /// Get efficiency of mesh distribution: In an ideal distribution
6066 /// without halo overhead, each processor would only hold its own
6067 /// elements. Efficieny per processor = (number of non-halo elements)/
6068 /// (total number of elements).
6069 //========================================================================
6071  double& max_efficiency,
6072  double& min_efficiency)
6073 {
6074  // Storage for number of processors and current processor
6075  int n_proc=Comm_pt->nproc();
6076  int my_rank=Comm_pt->my_rank();
6077 
6078  // Create vector to hold number of elements and halo elements
6079  Vector<int> nhalo_elements(n_proc);
6080  Vector<int> n_elements(n_proc);
6081 
6082  // Count total number of halo elements
6083  unsigned count=0;
6084  for (int d=0;d<n_proc;d++)
6085  {
6087  count+=halo_elem_pt.size();
6088  }
6089 
6090  // Stick own number into appropriate entry
6091  int nhalo_element_local=count;
6092  int n_element_local=nelement();
6093 
6094  // Gather information on root processor: First argument group
6095  // specifies what is to be sent (one int from each procssor, indicating
6096  // the number of elements on it), the second group indicates where
6097  // the results are to be gathered (in rank order) on root processor.
6098  MPI_Gather(&nhalo_element_local,1,MPI_INT,
6099  &nhalo_elements[0],1, MPI_INT,
6100  0,Comm_pt->mpi_comm());
6101  MPI_Gather(&n_element_local,1,MPI_INT,
6102  &n_elements[0],1, MPI_INT,
6103  0,Comm_pt->mpi_comm());
6104 
6105  // Initialise stats
6106  av_efficiency=0.0;
6107  double max=-1.0;
6108  double min=1000000000.0;
6109 
6110  if (my_rank==0)
6111  {
6112  for (int i=0;i<n_proc;i++)
6113  {
6114  unsigned nel=n_elements[i];
6115  double eff=0.0;
6116  if (nel!=0)
6117  {
6118  eff=double(n_elements[i]-nhalo_elements[i])/double(nel);
6119  }
6120  av_efficiency+=eff;
6121  if (eff>max) max=eff;
6122  if (eff<min) min=eff;
6123 
6124  }
6125  av_efficiency/=double(n_proc);
6126  }
6127 
6128  // Now broadcast the result back out
6129  MPI_Bcast(&max,1,MPI_DOUBLE,0,Comm_pt->mpi_comm());
6130  MPI_Bcast(&min,1,MPI_DOUBLE,0,Comm_pt->mpi_comm());
6131  MPI_Bcast(&av_efficiency,1,MPI_DOUBLE,0,Comm_pt->mpi_comm());
6132 
6133  max_efficiency=max;
6134  min_efficiency=min;
6135 
6136 }
6137 
6138 
6139 
6140 //========================================================================
6141 /// Doc the mesh distribution
6142 //========================================================================
6144 {
6145  // Storage for current processor and number of processors
6146  int my_rank=Comm_pt->my_rank();
6147  int n_proc=Comm_pt->nproc();
6148 
6149  std::ostringstream filename;
6150  std::ostringstream filename2;
6151  std::ofstream some_file;
6152  std::ofstream some_file2;
6153 
6154  // Doc elements on this processor
6155  filename << doc_info.directory() << "/"<< doc_info.label()<< "elements_on_proc"
6156  << my_rank << "_" << doc_info.number() << ".dat";
6157  some_file.open(filename.str().c_str());
6158  this->output(some_file,5);
6159  some_file.close();
6160 
6161  // Doc non-halo elements on this processor
6162  filename.str("");
6163  filename << doc_info.directory() << "/"
6164  << doc_info.label()<< "non_halo_elements_on_proc"
6165  << my_rank << "_" << doc_info.number() << ".dat";
6166  some_file.open(filename.str().c_str());
6167 
6168  // Get to elements on processor
6169  unsigned nelem=this->nelement();
6170  for (unsigned e=0;e<nelem;e++)
6171  {
6172  //Get the element
6173  GeneralisedElement* el_pt = this->element_pt(e);
6174  //Only output if not a halo element
6175  if(!el_pt->is_halo())
6176  {
6177  FiniteElement* f_el_pt=dynamic_cast<FiniteElement*>(el_pt);
6178  //Indicate a generalised element
6179  if(f_el_pt==0)
6180  {
6181  some_file << "Generalised Element " << e << "\n";
6182  }
6183  else
6184  // output if non-halo and a finite element
6185  {
6186  f_el_pt->output(some_file,5);
6187  }
6188  }
6189  }
6190 
6191  some_file.close();
6192 
6193 
6194  // Doc halo elements on this processor
6195  filename.str("");
6196  filename << doc_info.directory() << "/"<< doc_info.label()
6197  << "halo_elements_on_proc"
6198  << my_rank << "_" << doc_info.number() << ".dat";
6199  some_file.open(filename.str().c_str());
6200  for (int domain=0; domain<n_proc; domain++)
6201  {
6202  filename2.str("");
6203  filename2 << doc_info.directory() << "/"<< doc_info.label()
6204  << "halo_elements_with_proc" << domain << "_on_proc"
6205  << my_rank << "_" << doc_info.number() << ".dat";
6206  some_file2.open(filename2.str().c_str());
6207 
6208  // Get vector of halo elements by copy operation
6209  Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(domain));
6210  unsigned nelem=halo_elem_pt.size();
6211  for (unsigned e=0;e<nelem;e++)
6212  {
6213  FiniteElement* f_el_pt = dynamic_cast<FiniteElement*>(halo_elem_pt[e]);
6214  if(f_el_pt!=0)
6215  {
6216 
6217 #ifdef PARANOID
6218  // Check if it's refineable and if so if it's a leaf.
6219  // If not something must have gone wrong.
6220  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(
6221  f_el_pt);
6222  if (ref_el_pt!=0)
6223  {
6224  if (!(ref_el_pt->tree_pt()->is_leaf()))
6225  {
6226  std::ostringstream error_message;
6227  error_message
6228  << "Haloed element is not a leaf element. This shouldn't happen"
6229  << std::endl;
6230  error_message
6231  << "Here are the nodal positions: " << std::endl;
6232  unsigned nnod=ref_el_pt->nnode();
6233  for (unsigned j=0;j<nnod;j++)
6234  {
6235  Node* nod_pt=ref_el_pt->node_pt(j);
6236  unsigned n_dim=nod_pt->ndim();
6237  for (unsigned i=0;i<n_dim;i++)
6238  {
6239  error_message << nod_pt->x(i) << " ";
6240  }
6241  error_message << "\n";
6242  throw OomphLibError(error_message.str(),
6243  OOMPH_CURRENT_FUNCTION,
6244  OOMPH_EXCEPTION_LOCATION);
6245  }
6246  }
6247  }
6248 #endif
6249 
6250  f_el_pt->output(some_file,5);
6251  f_el_pt->output(some_file2,5);
6252  }
6253  //Indicate a generalised element
6254  else
6255  {
6256  some_file << "Generalised Element " << e << "\n";
6257  some_file2 << "Generalised Element " << e << "\n";
6258  }
6259  }
6260  some_file2.close();
6261  }
6262  some_file.close();
6263 
6264  // Doc haloed elements on this processor
6265  filename.str("");
6266  filename << doc_info.directory() << "/"<< doc_info.label()
6267  << "haloed_elements_on_proc"
6268  << my_rank << "_" << doc_info.number() << ".dat";
6269  some_file.open(filename.str().c_str());
6270  for (int domain=0; domain<n_proc; domain++)
6271  {
6272  filename2.str("");
6273  filename2 << doc_info.directory() << "/"<< doc_info.label()
6274  << "haloed_elements_with_proc" << domain << "_on_proc"
6275  << my_rank << "_" << doc_info.number() << ".dat";
6276  some_file2.open(filename2.str().c_str());
6277 
6278  // Get vector of haloed elements by copy operation
6279  Vector<GeneralisedElement*> haloed_elem_pt(this->haloed_element_pt(domain));
6280  unsigned nelem=haloed_elem_pt.size();
6281  for (unsigned e=0;e<nelem;e++)
6282  {
6283  //Is it a finite element
6284  FiniteElement *finite_el_pt = dynamic_cast<FiniteElement*>
6285  (haloed_elem_pt[e]);
6286  if(finite_el_pt!=0)
6287  {
6288 
6289 #ifdef PARANOID
6290  // Check if it's refineable and if so if it's a leaf.
6291  // If not something must have gone wrong.
6292  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(
6293  finite_el_pt);
6294  if (ref_el_pt!=0)
6295  {
6296  if (!(ref_el_pt->tree_pt()->is_leaf()))
6297  {
6298  std::ostringstream error_message;
6299  error_message
6300  << "Haloed element is not a leaf element. This shouldn't happen"
6301  << std::endl;
6302  error_message
6303  << "Here are the nodal positions: " << std::endl;
6304  unsigned nnod=ref_el_pt->nnode();
6305  for (unsigned j=0;j<nnod;j++)
6306  {
6307  Node* nod_pt=ref_el_pt->node_pt(j);
6308  unsigned n_dim=nod_pt->ndim();
6309  for (unsigned i=0;i<n_dim;i++)
6310  {
6311  error_message << nod_pt->x(i) << " ";
6312  }
6313  error_message << "\n";
6314  throw OomphLibError(error_message.str(),
6315  OOMPH_CURRENT_FUNCTION,
6316  OOMPH_EXCEPTION_LOCATION);
6317  }
6318  }
6319  }
6320 #endif
6321 
6322  finite_el_pt->output(some_file,5);
6323  finite_el_pt->output(some_file2,5);
6324  }
6325  //Indicate a generalised element
6326  else
6327  {
6328  some_file << "Generalised Element " << e << "\n";
6329  some_file2 << "Generalised Element " << e << "\n";
6330  }
6331  }
6332  some_file2.close();
6333  }
6334  some_file.close();
6335 
6336 
6337  // Doc non-halo nodes on this processor
6338  filename.str("");
6339  filename << doc_info.directory() << "/"<< doc_info.label()
6340  << "non_halo_nodes_on_proc"
6341  << my_rank << "_" << doc_info.number() << ".dat";
6342  some_file.open(filename.str().c_str());
6343  unsigned nnod=this->nnode();
6344  for (unsigned j=0;j<nnod;j++)
6345  {
6346  Node* nod_pt=this->node_pt(j);
6347  if (!nod_pt->is_halo())
6348  {
6349  unsigned ndim=nod_pt->ndim();
6350  for (unsigned i=0;i<ndim;i++)
6351  {
6352  some_file << nod_pt->x(i) << " " ;
6353  }
6354  some_file << nod_pt->nvalue() << " "
6355  << nod_pt->non_halo_proc_ID() << " "
6356  << nod_pt->is_halo() << " "
6357  << nod_pt->hang_code() << std::endl;
6358  }
6359  }
6360  some_file.close();
6361 
6362 
6363  // Doc nodes on this processor
6364  filename.str("");
6365  filename << doc_info.directory() << "/"<< doc_info.label()
6366  << "nodes_on_proc"
6367  << my_rank << "_" << doc_info.number() << ".dat";
6368  some_file.open(filename.str().c_str());
6369  for (unsigned j=0;j<nnod;j++)
6370  {
6371  Node* nod_pt=this->node_pt(j);
6372  unsigned ndim=nod_pt->ndim();
6373  for (unsigned i=0;i<ndim;i++)
6374  {
6375  some_file << nod_pt->x(i) << " " ;
6376  }
6377  some_file << nod_pt->nvalue() << " "
6378  << nod_pt->non_halo_proc_ID() << " "
6379  << nod_pt->is_halo() << " "
6380  << nod_pt->hang_code() << std::endl;
6381  }
6382  some_file.close();
6383 
6384  // Doc solid nodes on this processor
6385  filename.str("");
6386  filename << doc_info.directory() << "/"
6387  << doc_info.label()<< "solid_nodes_on_proc"
6388  << my_rank << "_" << doc_info.number() << ".dat";
6389  some_file.open(filename.str().c_str());
6390  unsigned nsnod=this->nnode();
6391  for (unsigned j=0;j<nsnod;j++)
6392  {
6393  Node* nod_pt=this->node_pt(j);
6394  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
6395  if (solid_nod_pt!=0)
6396  {
6397  unsigned ndim=solid_nod_pt->ndim();
6398  for (unsigned i=0;i<ndim;i++)
6399  {
6400  some_file << nod_pt->x(i) << " " ;
6401  }
6402  some_file << nod_pt->nvalue() << " "
6403  << nod_pt->non_halo_proc_ID() << " "
6404  << nod_pt->is_halo() << " "
6405  << nod_pt->hang_code() << "\n";
6406  }
6407  }
6408  some_file.close();
6409 
6410 
6411  // Doc halo nodes on this processor
6412  filename.str("");
6413  filename << doc_info.directory() << "/"<< doc_info.label()
6414  << "halo_nodes_on_proc"
6415  << my_rank << "_" << doc_info.number() << ".dat";
6416  some_file.open(filename.str().c_str());
6417  for (int domain=0; domain<n_proc; domain++)
6418  {
6419  filename2.str("");
6420  filename2 << doc_info.directory() << "/"<< doc_info.label()
6421  << "halo_nodes_with_proc" << domain << "_on_proc"
6422  << my_rank << "_" << doc_info.number() << ".dat";
6423  some_file2.open(filename2.str().c_str());
6424  unsigned nnod=this->nhalo_node(domain);
6425  for (unsigned j=0;j<nnod;j++)
6426  {
6427  Node* nod_pt=this->halo_node_pt(domain,j);
6428  unsigned ndim=nod_pt->ndim();
6429  for (unsigned i=0;i<ndim;i++)
6430  {
6431  some_file << nod_pt->x(i) << " " ;
6432  some_file2 << nod_pt->x(i) << " " ;
6433  }
6434  some_file << nod_pt->nvalue() << " "
6435  << nod_pt->non_halo_proc_ID() << " "
6436  << nod_pt->is_halo() << " "
6437  << nod_pt->hang_code() << "\n";
6438  some_file2 << nod_pt->nvalue() << " "
6439  << nod_pt->non_halo_proc_ID() << " "
6440  << nod_pt->is_halo() << " "
6441  << nod_pt->hang_code() << "\n";
6442  }
6443  some_file2.close();
6444  }
6445  some_file.close();
6446 
6447 
6448 
6449  // Doc haloed nodes on this processor
6450  filename.str("");
6451  filename << doc_info.directory() << "/"<< doc_info.label()
6452  << "haloed_nodes_on_proc"
6453  << my_rank << "_" << doc_info.number() << ".dat";
6454  some_file.open(filename.str().c_str());
6455  for (int domain=0; domain<n_proc; domain++)
6456  {
6457  filename2.str("");
6458  filename2 << doc_info.directory() << "/"<< doc_info.label()
6459  << "haloed_nodes_with_proc" << domain << "_on_proc"
6460  << my_rank << "_" << doc_info.number() << ".dat";
6461  some_file2.open(filename2.str().c_str());
6462  unsigned nnod=this->nhaloed_node(domain);
6463  for (unsigned j=0;j<nnod;j++)
6464  {
6465  Node* nod_pt=this->haloed_node_pt(domain,j);
6466  unsigned ndim=nod_pt->ndim();
6467  for (unsigned i=0;i<ndim;i++)
6468  {
6469  some_file << nod_pt->x(i) << " " ;