mesh.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 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 //Header file for general mesh classes
31 
32 //Include guards to prevent multiple inclusion of the header
33 #ifndef OOMPH_GENERIC_MESH_HEADER
34 #define OOMPH_GENERIC_MESH_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38  #include <oomph-lib-config.h>
39 #endif
40 
41 #ifdef OOMPH_HAS_MPI
42 #include "mpi.h"
43 #endif
44 
45 #include<float.h>
46 #include <list>
47 #include <typeinfo>
48 #include<string>
49 
50 //oomph-lib headers
51 #include "Vector.h"
52 #include "nodes.h"
53 #include "elements.h"
54 #include "timesteppers.h"
56 #include "matrices.h"
57 #include "refineable_elements.h"
58 
59 namespace oomph
60 {
61 
62 
63 
64 
65 //=================================================================
66 /// \short A general mesh class.
67 ///
68 /// The main components of a Mesh are:
69 /// - pointers to its Nodes
70 /// - pointers to its Elements
71 /// - pointers to its boundary Nodes
72 ///
73 //=================================================================
74 class Mesh
75 {
76 
77  public:
78 
79  /// Problem is a friend
80  friend class Problem;
81 
82 
83  /// \short Default Steady Timestepper, to be used in default arguments
84  /// to Mesh constructors
86 
87 
88  protected:
89 
90  /// \short Vector of Vector of pointers to nodes on the boundaries:
91  /// Boundary_node_pt(b,n). Note that this is private to force
92  /// the use of the add_boundary_node() function, which ensures
93  /// that the reverse look-up schemes for the nodes are set up.
95 
96  /// Flag to indicate that the lookup schemes for elements that are adjacent
97  /// to the boundaries has been set up.
99 
100  /// \short Vector of Vector of pointers to elements on the boundaries:
101  /// Boundary_element_pt(b,e)
103 
104  /// \short For the e-th finite element on boundary b, this is the index of
105  /// the face that lies along that boundary
107 
108 #ifdef OOMPH_HAS_MPI
109 
110  /// Map of vectors holding the pointers to the root halo elements
111  std::map<unsigned, Vector<GeneralisedElement*> > Root_halo_element_pt;
112 
113  /// Map of vectors holding the pointers to the root haloed elements
114  std::map<unsigned, Vector<GeneralisedElement*> > Root_haloed_element_pt;
115 
116  /// Map of vectors holding the pointers to the halo nodes
117  std::map<unsigned, Vector<Node*> > Halo_node_pt;
118 
119  /// Map of vectors holding the pointers to the haloed nodes
120  std::map<unsigned, Vector<Node*> > Haloed_node_pt;
121 
122  /// Map of vectors holding the pointers to the shared nodes.
123  /// These are all the nodes that are on two "neighbouring" processes
124  /// (the halo(ed) lookup scheme depends upon which processor is in charge
125  /// - a node which is on 3 processors, for example, will not feature in
126  /// the halo(ed) lookup scheme between the two lowest-numbered processors)
127  std::map<unsigned, Vector<Node*> > Shared_node_pt;
128 
129  /// Pointer to communicator -- set to NULL if mesh is not distributed
131 
132  /// External halo(ed) elements are created as and when they are needed
133  /// to act as source elements for the particular process's mesh.
134  /// The storage is wiped and rebuilt every time the mesh is refined.
135 
136  /// Map of vectors holding the pointers to the external halo elements
137  std::map<unsigned, Vector<GeneralisedElement*> > External_halo_element_pt;
138 
139  /// Map of vectors holding the pointers to the external haloed elements
140  std::map<unsigned, Vector<GeneralisedElement*> > External_haloed_element_pt;
141 
142 
143  // External halo(ed) nodes are on the external halo(ed) elements
144 
145  /// Map of vectors holding the pointers to the external halo nodes
146  std::map<unsigned, Vector<Node*> > External_halo_node_pt;
147 
148  /// Map of vectors holding the pointers to the external haloed nodes
149  std::map<unsigned, Vector<Node*> > External_haloed_node_pt;
150 
151  /// bool to indicate whether to keep all elements in a mesh as halos or not
153 
154  /// Set this to true to suppress resizing of halo nodes (at your own risk!)
156 
157  /// Setup shared node scheme
159 
160 #endif
161 
162  /// \short Assign the global equation numbers in the Data stored at the nodes
163  /// and also internal element Data. Also, build (via push_back) the
164  /// Vector of pointers to the dofs (variables).
165  unsigned long
167 
168  /// \short Function to describe the dofs of the Mesh. The ostream
169  /// specifies the output stream to which the description
170  /// is written; the string stores the currently
171  /// assembled output that is ultimately written to the
172  /// output stream by Data::describe_dofs(...); it is typically
173  /// built up incrementally as we descend through the
174  /// call hierarchy of this function when called from
175  /// Problem::describe_dofs(...)
176  void describe_dofs(std::ostream& out,const std::string& current_string) const;
177 
178  /// \short Function to describe the local dofs of the elements. The ostream
179  /// specifies the output stream to which the description
180  /// is written; the string stores the currently
181  /// assembled output that is ultimately written to the
182  /// output stream by Data::describe_dofs(...); it is typically
183  /// built up incrementally as we descend through the
184  /// call hierarchy of this function when called from
185  /// Problem::describe_dofs(...)
186  void describe_local_dofs(std::ostream& out,
187  const std::string& current_string) const;
188 
189  /// \short Assign the local equation numbers in all elements
190  /// If the boolean argument is true then also store pointers to dofs
191  void assign_local_eqn_numbers(const bool &store_local_dof_pt);
192 
193  /// Vector of pointers to nodes
195 
196  /// Vector of pointers to generalised elements
198 
199  /// \short Vector of boolean data that indicates whether the boundary
200  /// coordinates have been set for the boundary
201  std::vector<bool> Boundary_coordinate_exists;
202 
203  /// \short A function that upgrades an ordinary node to a boundary node
204  /// We shouldn't ever really use this, but it does make life that
205  /// bit easier for the lazy mesh writer. The pointer to the node is
206  /// replaced by a pointer to the new boundary node in all element look-up
207  /// schemes and in the mesh's Node_pt vector. The new node is also
208  /// addressed by node_pt on return from the function.
211 
213 
214 
215 public:
216 
217 
218 #ifdef OOMPH_HAS_MPI
219 
220 
221  /// \short Helper function that resizes halo nodes to the same
222  /// size as their non-halo counterparts if required. (A discrepancy
223  /// can arise if a FaceElement that introduces additional unknowns
224  /// are attached to a bulk element that shares a node with a haloed element.
225  /// In that case the joint node between haloed and non-haloed element
226  /// is resized on that processor but not on the one that holds the
227  /// halo counterpart (because no FaceElement is attached to the halo
228  /// element)
229  void resize_halo_nodes();
230 
231 #endif
232 
233 
234  /// \short Typedef for function pointer to function that computes
235  /// steady exact solution
237  const Vector<double>& x, Vector<double>& soln);
238 
239  /// \short Typedef for function pointer to function that computes unsteady
240  /// exact solution
242  const double& time, const Vector<double>& x, Vector<double>& soln);
243 
244  /// \short Boolean used to control warning about empty mesh level
245  /// timestepper function
247 
248  /// \short Default constructor
250  {
251  // Lookup scheme hasn't been setup yet
253 #ifdef OOMPH_HAS_MPI
254  // Set defaults for distributed meshes
255 
256  // Mesh hasn't been distributed: Null out pointer to communicator
257  Comm_pt=0;
258  // Don't keep all objects as halos
260  // Don't output halo elements
261  Output_halo_elements=false;
262  // Don't suppress automatic resizing of halo nodes
264 #endif
265  }
266 
267  /// \short Constructor builds combined mesh from the meshes specified.
268  /// Note: This simply merges the meshes' elements and nodes (ignoring
269  /// duplicates; no boundary information etc. is created).
270  Mesh(const Vector<Mesh*>& sub_mesh_pt)
271  {
272 #ifdef OOMPH_HAS_MPI
273  // Mesh hasn't been distributed: Null out pointer to communicator
274  Comm_pt=0;
275 #endif
276  // Now merge the meshes
277  merge_meshes(sub_mesh_pt);
278  }
279 
280  /// \short Merge meshes.
281  /// Note: This simply merges the meshes' elements and nodes (ignoring
282  /// duplicates; no boundary information etc. is created).
283  void merge_meshes(const Vector<Mesh*>& sub_mesh_pt);
284 
285  /// \short Interface for function that is used to setup the boundary
286  /// information (Empty virtual function -- implement this for specific
287  /// Mesh classes)
288  virtual void setup_boundary_element_info() { }
289 
290  /// \short Setup lookup schemes which establish whic elements are located
291  /// next to mesh's boundaries. Doc in outfile (if it's open).
292  /// (Empty virtual function -- implement this for specific
293  /// Mesh classes)
294  virtual void setup_boundary_element_info(std::ostream &outfile) {}
295 
296  /// Virtual function to perform the reset boundary elements info rutines
298  Vector<unsigned> &ntmp_boundary_elements,
299  Vector<Vector<unsigned> > &ntmp_boundary_elements_in_region,
300  Vector<FiniteElement*> &deleted_elements)
301  {
302  std::ostringstream error_stream;
303  error_stream << "Empty default reset boundary element info function"
304  << "called.\n";
305  error_stream << "This should be overloaded in a specific "
306  << "TriangleMeshBase\n";
307  throw OomphLibError(error_stream.str(),
308  "Mesh::reset_boundary_element_info()",
309  OOMPH_EXCEPTION_LOCATION);
310  }
311 
312  /// \short Output boundary coordinates on boundary b -- template argument
313  /// specifies the bulk element type (needed to create FaceElement
314  /// of appropriate type on mesh boundary).
315  template<class BULK_ELEMENT>
316  void doc_boundary_coordinates(const unsigned& b, std::ofstream& the_file)
317  {
318  if (nelement()==0) return;
320  {
321  oomph_info << "No boundary coordinates were set up for boundary "
322  << b << std::endl;
323  return;
324  }
325 
326  // Get spatial dimension
327  unsigned dim=finite_element_pt(0)->node_pt(0)->ndim();
328 
329  // Loop over all elements on boundaries
330  unsigned nel=this->nboundary_element(b);
331 
332  // Loop over the bulk elements adjacent to boundary b
333  for(unsigned e=0;e<nel;e++)
334  {
335  // Get pointer to the bulk element that is adjacent to boundary b
336  FiniteElement* bulk_elem_pt = this->boundary_element_pt(b,e);
337 
338  //Find the index of the face of element e along boundary b
339  int face_index = this->face_index_at_boundary(b,e);
340 
341  // Create new face element
343  new DummyFaceElement<BULK_ELEMENT>(bulk_elem_pt,face_index);
344 
345  // Specify boundary id in bulk mesh (needed to extract
346  // boundary coordinate)
348 
349  // Doc boundary coordinate
350  Vector<double> s(dim-1);
351  Vector<double> zeta(dim-1);
352  Vector<double> x(dim);
353  unsigned n_plot=5;
354  the_file << el_pt->tecplot_zone_string(n_plot);
355 
356  // Loop over plot points
357  unsigned num_plot_points=el_pt->nplot_points(n_plot);
358  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
359  {
360  // Get local coordinates of plot point
361  el_pt->get_s_plot(iplot,n_plot,s);
362  el_pt->interpolated_zeta(s,zeta);
363  el_pt->interpolated_x(s,x);
364  for (unsigned i=0;i<dim;i++)
365  {
366  the_file << x[i] << " ";
367  }
368  for (unsigned i=0;i<(dim-1);i++)
369  {
370  the_file << zeta[i] << " ";
371  }
372 
373  the_file << std::endl;
374  }
375  el_pt->write_tecplot_zone_footer(the_file,n_plot);
376 
377  // Cleanup
378  delete el_pt;
379  }
380  }
381 
382 
383  /// \short Scale all nodal coordinates by given factor. Virtual
384  /// so it can be overloaded in SolidMesh class where it also
385  /// re-assigns the Lagrangian coordinates.
386  virtual void scale_mesh(const double& factor)
387  {
388  unsigned nnod=this->nnode();
389  unsigned dim=this->node_pt(0)->ndim();
390  for (unsigned j=0;j<nnod;j++)
391  {
392  Node* nod_pt=this->node_pt(j);
393  for (unsigned i=0;i<dim;i++)
394  {
395  nod_pt->x(i)*=factor;
396  }
397  }
398  }
399 
400 
401  /// Broken copy constructor
402  Mesh(const Mesh& dummy)
403  {
404  BrokenCopy::broken_copy("Mesh");
405  }
406 
407  /// Broken assignment operator
408  void operator=(const Mesh&)
409  {
411  }
412 
413  /// Virtual Destructor to clean up all memory
414  virtual ~Mesh();
415 
416 
417 /// \short Flush storage for elements and nodes by emptying the
418 /// vectors that store the pointers to them. This is
419 /// useful if a particular mesh is only built to generate
420 /// a small part of a bigger mesh. Once the elements and
421 /// nodes have been created, they are typically copied
422 /// into the new mesh and the auxiliary mesh can be
423 /// deleted. However, if we simply call the destructor
424 /// of the auxiliary mesh, it will also wipe out
425 /// the nodes and elements, because it still "thinks"
426 /// it's in charge of these...
428  {
431  }
432 
433 /// \short Flush storage for elements (only) by emptying the
434 /// vectors that store the pointers to them. This is
435 /// useful if a particular mesh is only built to generate
436 /// a small part of a bigger mesh. Once the elements and
437 /// nodes have been created, they are typically copied
438 /// into the new mesh and the auxiliary mesh can be
439 /// deleted. However, if we simply call the destructor
440 /// of the auxiliary mesh, it will also wipe out
441 /// the nodes and elements, because it still "thinks"
442 /// it's in charge of these...
444  {
445  Element_pt.clear();
446  }
447 
448 /// \short Flush storage for nodes (only) by emptying the
449 /// vectors that store the pointers to them.
451  {
452  Node_pt.clear();
453  }
454 
455  /// Return pointer to global node n
456  Node* &node_pt(const unsigned long &n) {return Node_pt[n];}
457 
458  /// Return pointer to global node n (const version)
459  Node* node_pt(const unsigned long &n) const {return Node_pt[n];}
460 
461  /// Return pointer to element e
462  GeneralisedElement* &element_pt(const unsigned long &e)
463  {return Element_pt[e];}
464 
465  /// Return pointer to element e (const version)
466  GeneralisedElement* element_pt(const unsigned long &e) const
467  {return Element_pt[e];}
468 
469  /// Return reference to the Vector of elements
471 
472  /// Return reference to the Vector of elements
474 
475  /// \short Upcast (downcast?) to FiniteElement
476  /// (needed to access FiniteElement member functions).
477  FiniteElement* finite_element_pt(const unsigned &e) const
478  {
479 #ifdef PARANOID
480  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
481  if (el_pt==0)
482  {
483  //Error
484  throw OomphLibError("Failed cast to FiniteElement* ",
485  OOMPH_CURRENT_FUNCTION,
486  OOMPH_EXCEPTION_LOCATION);
487  //Dummy return to keep intel compiler happy
488  return el_pt;
489  }
490  return el_pt;
491 #else
492  return dynamic_cast<FiniteElement*>(Element_pt[e]);
493 #endif
494  }
495 
496  /// Return pointer to node n on boundary b
497  Node* &boundary_node_pt(const unsigned &b, const unsigned &n)
498  {return Boundary_node_pt[b][n];}
499 
500  /// Return pointer to node n on boundary b
501  Node* boundary_node_pt(const unsigned &b, const unsigned &n) const
502  {return Boundary_node_pt[b][n];}
503 
504  /// Set the number of boundaries in the mesh
505  void set_nboundary(const unsigned &nbound)
506  {
507  Boundary_node_pt.resize(nbound);
508  Boundary_coordinate_exists.resize(nbound,false);
509  }
510 
511  ///\short Clear all pointers to boundary nodes
512  void remove_boundary_nodes();
513 
514  /// \short Remove all information about nodes stored on the b-th
515  /// boundary of the mesh
516  void remove_boundary_nodes(const unsigned &b);
517 
518  //\ short Remove a node from the boundary b
519  void remove_boundary_node(const unsigned &b, Node* const &node_pt);
520 
521  /// Add a (pointer to) a node to the b-th boundary
522  void add_boundary_node(const unsigned &b, Node* const &node_pt);
523 
524  /// Replace existing boundary node lookup schemes with new schemes created
525  /// using the boundary data stored in the nodes.
527  {
528  // Clear existing boundary data
529  Boundary_node_pt.clear();
530 
531  // Loop over nodes adding them to the appropriate boundary lookup schemes
532  // in the mesh.
533  const unsigned n_node = nnode();
534  for(unsigned nd=0; nd<n_node; nd++)
535  {
536  Node* nd_pt = node_pt(nd);
537 
538  if(nd_pt->is_on_boundary())
539  {
540  // Get set of boundaries that the node is on
541  std::set<unsigned>* boundaries_pt;
542  nd_pt->get_boundaries_pt(boundaries_pt);
543 
544  // If needed then add more boundaries to this mesh
545  unsigned max_boundary_n = 1 +
546  *std::max_element(boundaries_pt->begin(), boundaries_pt->end());
547  if(max_boundary_n > nboundary())
548  {
549  set_nboundary(max_boundary_n);
550  }
551 
552  // Add node pointer to the appropriate Boundary_node_pt vectors
553  std::set<unsigned>::const_iterator it;
554  for(it=boundaries_pt->begin(); it!=boundaries_pt->end(); it++)
555  {
556  Boundary_node_pt[*it].push_back(nd_pt);
557  }
558  }
559  }
560  }
561 
562  /// Indicate whether the i-th boundary has an intrinsic coordinate.
563  // By default, if the array Boundary_coordinate has not been resized,
564  // return false
565  bool boundary_coordinate_exists(const unsigned &i) const
566  {
567  if(Boundary_coordinate_exists.empty()) {return false;}
568  //ALH: This bounds-checking code needs to remain, because
569  //Boundary_coordinate_exists is
570  //an stl vector not our overloaded Vector class
571 #ifdef RANGE_CHECKING
572  if (i>=Boundary_coordinate_exists.size())
573  {
574  std::ostringstream error_message;
575  error_message << "Range Error: " << i << " is not in the range (0,"
576  << Boundary_coordinate_exists.size()-1 << std::endl;
577 
578  throw OomphLibError(error_message.str(),
579  OOMPH_CURRENT_FUNCTION,
580  OOMPH_EXCEPTION_LOCATION);
581  }
582 #endif
584  }
585 
586  /// Return number of elements in the mesh
587  unsigned long nelement() const {return Element_pt.size();}
588 
589  /// Return number of nodes in the mesh
590  unsigned long nnode() const {return Node_pt.size();}
591 
592  /// Return number of dof types in mesh
593  unsigned ndof_types() const;
594 
595  /// Return number of elemental dimension in mesh
596  unsigned elemental_dimension() const;
597 
598  /// Return number of nodal dimension in mesh
599  unsigned nodal_dimension() const;
600 
601  /// Add a (pointer to a) node to the mesh
602  void add_node_pt(Node* const &node_pt) {Node_pt.push_back(node_pt);}
603 
604  /// Add a (pointer to) an element to the mesh
606  {Element_pt.push_back(element_pt);}
607 
608  /// \short Update nodal positions in response to changes in the domain shape.
609  /// Uses the FiniteElement::get_x(...) function for FiniteElements
610  /// and doesn't do anything for other element types.
611  /// If a MacroElement pointer has been set for a FiniteElement,
612  /// the MacroElement representation is used to update the
613  /// nodal positions; if not get_x(...) uses the FE interpolation
614  /// and thus leaves the nodal positions unchanged.
615  /// Virtual, so it can be overloaded by specific meshes,
616  /// such as AlgebraicMeshes or SpineMeshes.
617  /// Generally, this function updates the position of all nodes
618  /// in response to changes in the boundary position.
619  /// However, we ignore all SolidNodes since their
620  /// position is computed as part of the solution -- unless
621  /// the bool flag is set to true. Such calls are typically made
622  /// when the initial mesh is created and/or after a mesh has been
623  /// refined repeatedly before the start of the computation.
624  virtual void node_update(const bool& update_all_solid_nodes=false);
625 
626  /// \short Re-order nodes in the order in which they appear in elements --
627  /// can be overloaded for more efficient re-ordering
628  virtual void reorder_nodes(const bool& use_old_ordering=true);
629 
630  /// \short Get a reordering of the nodes in the order in which they
631  /// appear in elements -- can be overloaded for more efficient
632  /// re-ordering
633  virtual void get_node_reordering(Vector<Node*> &reordering,
634  const bool& use_old_ordering=true) const;
635 
636  /// \short Constuct a Mesh of FACE_ELEMENTs along the b-th boundary
637  /// of the mesh (which contains elements of type BULK_ELEMENT)
638  template<class BULK_ELEMENT, template<class> class FACE_ELEMENT>
639  void build_face_mesh(const unsigned &b, Mesh* const &face_mesh_pt)
640  {
641  //Find the number of nodes on the boundary
642  unsigned nbound_node = nboundary_node(b);
643  //Loop over the boundary nodes and add them to face mesh node pointer
644  for(unsigned n=0;n<nbound_node;n++)
645  {face_mesh_pt->add_node_pt(boundary_node_pt(b,n));}
646 
647  //Find the number of elements next to the boundary
648  unsigned nbound_element = nboundary_element(b);
649  //Loop over the elements adjacent to boundary b
650  for(unsigned e=0;e<nbound_element;e++)
651  {
652  //Create the FaceElement
653  FACE_ELEMENT<BULK_ELEMENT>* face_element_pt =
654  new FACE_ELEMENT<BULK_ELEMENT>(boundary_element_pt(b,e),
656 
657  //Add the face element to the face mesh
658  face_mesh_pt->add_element_pt(face_element_pt);
659  }
660 
661 #ifdef OOMPH_HAS_MPI
662  // If the bulk mesh has been distributed then the face mesh is too
663  if (this->is_mesh_distributed())
664  {
665  face_mesh_pt->set_communicator_pt(this->communicator_pt());
666  }
667 #endif
668  }
669 
670  /// \short Self-test: Check elements and nodes. Return 0 for OK
671  unsigned self_test();
672 
673 
674  /// \short Determine max and min area for all FiniteElements in the mesh
675  /// (non-FiniteElements are ignored)
676  void max_and_min_element_size(double& max_size, double& min_size)
677  {
678  max_size=0.0;
679  min_size=DBL_MAX;
680  unsigned nel=nelement();
681  for (unsigned e=0;e<nel;e++)
682  {
683  max_size=std::max(max_size,
684  finite_element_pt(e)->size());
685  min_size=std::min(min_size,
686  finite_element_pt(e)->size()); }
687  }
688 
689 
690  /// \short Determine the sum of all "sizes" of the FiniteElements in the mesh
691  /// (non-FiniteElements are ignored). This gives the length/area/volume
692  /// occupied by the mesh
693  double total_size()
694  {
695  double size=0.0;
696  unsigned nel=nelement();
697  for (unsigned e=0;e<nel;e++)
698  {
700  if (fe_pt!=0)
701  {
702  size+=fe_pt->size();
703  }
704  }
705  return size;
706  }
707 
708 
709  /// \short Check for inverted elements and report outcome
710  /// in boolean variable. This visits all elements at their
711  /// integration points and checks if the Jacobian of the
712  /// mapping between local and global coordinates is positive --
713  /// using the same test that would be carried out (but only in PARANOID
714  /// mode) during the assembly of the elements' Jacobian matrices.
715  /// Inverted elements are output in inverted_element_file (if the
716  /// stream is open).
717  void check_inverted_elements(bool& mesh_has_inverted_elements,
718  std::ofstream& inverted_element_file);
719 
720 
721 /// \short Check for inverted elements and report outcome
722  /// in boolean variable. This visits all elements at their
723  /// integration points and checks if the Jacobian of the
724  /// mapping between local and global coordinates is positive --
725  /// using the same test that would be carried out (but only in PARANOID
726  /// mode) during the assembly of the elements' Jacobian matrices.
727  void check_inverted_elements(bool& mesh_has_inverted_elements)
728  {
729  std::ofstream inverted_element_file;
730  check_inverted_elements(mesh_has_inverted_elements,
731  inverted_element_file);
732  }
733 
734 
735  /// \short Check for repeated nodes within a given spatial tolerance.
736  /// Return (0/1) for (pass/fail).
737  unsigned check_for_repeated_nodes(const double& epsilon=1.0e-12)
738  {
739  oomph_info <<"\n\nStarting check for repeated nodes...";
740  bool failed=false;
741  unsigned nnod=nnode();
742  for (unsigned j=0;j<nnod;j++)
743  {
744  Node* nod1_pt=this->node_pt(j);
745  unsigned dim=nod1_pt->ndim();
746  for (unsigned k=j+1;k<nnod;k++)
747  {
748  Node* nod2_pt=this->node_pt(k);
749  double dist=0.0;
750  for (unsigned i=0;i<dim;i++)
751  {
752  dist+=pow((nod1_pt->x(i)-nod2_pt->x(i)),2);
753  }
754  dist=sqrt(dist);
755  if (dist<epsilon)
756  {
757  oomph_info << "\n\nRepeated node!" << std::endl;
758  oomph_info << "Distance between nodes " << j << " and " << k
759  << std::endl;
760  oomph_info << "is " << dist << " which is less than the"
761  << std::endl;
762  oomph_info << "permitted distance of " << epsilon
763  << std::endl << std::endl;
764  oomph_info << "The offending nodes are located at: " << std::endl;
765  for (unsigned i=0;i<dim;i++)
766  {
767  oomph_info << nod1_pt->x(i) << " ";
768  }
769  if (nod1_pt->is_a_copy()||
770  nod2_pt->is_a_copy())
771  {
772  oomph_info
773  << "\n\n[NOTE: message issued as diagonistic rather than an error\n"
774  << " because at least one of the nodes is a copy; you may still\n"
775  << " want to check this out. BACKGROUND: Copied nodes share the same Data but\n"
776  << " will, in general, have different spatial positions (e.g. when used\n"
777  << " as periodic nodes); however there are cases when they are located\n"
778  << " at the same spatial position (e.g. in oomph-lib's annular mesh which\n"
779  << " is a rolled-around version of the rectangular quadmesh). In such cases,\n"
780  << " the nodes could have been deleted and completely replaced by \n"
781  << " pointers to existing nodes, but may have been left there for convenience\n"
782  << " or out of laziness...]\n";
783  }
784  else
785  {
786  failed=true;
787  }
788  oomph_info << std::endl << std::endl;
789 
790  }
791  }
792  }
793  if (failed) return 1;
794 
795  // If we made it to here, we must have passed the test.
796  oomph_info <<"...done: Test passed!" << std::endl << std::endl;
797  return 0;
798  }
799 
800  /// \short Prune nodes. Nodes that have been marked as obsolete are removed
801  /// from the mesh (and its boundary-node scheme). Returns vector
802  /// of pointers to deleted nodes.
804 
805  /// Return number of boundaries
806  unsigned nboundary() const {return Boundary_node_pt.size();}
807 
808  /// Return number of nodes on a particular boundary
809  unsigned long nboundary_node(const unsigned &ibound) const
810  {return Boundary_node_pt[ibound].size();}
811 
812 
813  /// Return pointer to e-th finite element on boundary b
814  FiniteElement* boundary_element_pt(const unsigned &b, const unsigned &e) const
815  {
816 #ifdef PARANOID
818  {
819  throw OomphLibError(
820  "Lookup scheme for elements next to boundary hasn't been set up yet!\n",
821  OOMPH_CURRENT_FUNCTION,
822  OOMPH_EXCEPTION_LOCATION);
823  }
824 #endif
825  return Boundary_element_pt[b][e];
826  }
827 
828 
829  /// \short Find a node not on any boundary in mesh_pt (useful for pinning
830  /// a single node in a purely Neumann problem so that it is fully
831  /// determined).
833  {
834  for(unsigned nd=0, nnd=nnode(); nd<nnd; nd++)
835  {
836  if( !(node_pt(nd)->is_on_boundary()))
837  {
838  return node_pt(nd);
839  }
840  }
841 
842  std::ostringstream error_msg;
843  error_msg << "No non-boundary nodes in the mesh.";
844  throw OomphLibError(error_msg.str(),
845  OOMPH_CURRENT_FUNCTION,
846  OOMPH_EXCEPTION_LOCATION);
847  // Never get here!
848  return 0;
849  }
850 
851  /// Return number of finite elements that are adjacent to boundary b
852  unsigned nboundary_element(const unsigned &b) const
853  {
854 #ifdef PARANOID
856  {
857  throw OomphLibError(
858  "Lookup scheme for elements next to boundary hasn't been set up yet!\n",
859  OOMPH_CURRENT_FUNCTION,
860  OOMPH_EXCEPTION_LOCATION);
861  }
862 #endif
863  return Boundary_element_pt[b].size();
864  }
865 
866 
867  /// \short For the e-th finite element on boundary b, return int to indicate
868  /// the face_index of the face adjacent to the boundary. This is consistent
869  /// with input required during the generation of FaceElements.
870  int face_index_at_boundary(const unsigned &b, const unsigned &e) const
871  {
872 #ifdef PARANOID
874  {
875  throw OomphLibError(
876  "Lookup scheme for elements next to boundary hasn't been set up yet!\n",
877  OOMPH_CURRENT_FUNCTION,
878  OOMPH_EXCEPTION_LOCATION);
879  }
880 #endif
881  return Face_index_at_boundary[b][e];
882  }
883 
884  /// Dump the data in the mesh into a file for restart
885  virtual void dump(std::ofstream &dump_file,
886  const bool& use_old_ordering=true) const;
887 
888  /// Dump the data in the mesh into a file for restart
889  void dump(const std::string &dump_file_name,
890  const bool& use_old_ordering=true) const
891  {
892  std::ofstream dump_stream(dump_file_name.c_str());
893 #ifdef PARANOID
894  if(!dump_stream.is_open())
895  {
896  std::string err = "Couldn't open file "+dump_file_name;
897  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
898  OOMPH_CURRENT_FUNCTION);
899  }
900 #endif
901  dump(dump_stream, use_old_ordering);
902  }
903 
904  /// \short Read solution from restart file
905  virtual void read(std::ifstream &restart_file);
906 
907 
908  /// \short Output in paraview format into specified file. Breaks up each
909  /// element into sub-elements for plotting purposes. We assume
910  /// that all elements are of the same type (fct will break
911  /// break (in paranoid mode) if paraview output fcts of the
912  /// elements are inconsistent).
913  void output_paraview(std::ofstream &file_out,
914  const unsigned &nplot) const;
915 
916  /// \short Output in paraview format into specified file. Breaks up each
917  /// element into sub-elements for plotting purposes. We assume
918  /// that all elements are of the same type (fct will break
919  /// break (in paranoid mode) if paraview output fcts of the
920  /// elements are inconsistent).
921  void output_fct_paraview(std::ofstream &file_out,
922  const unsigned &nplot,
923  FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const;
924 
925  /// Output for all elements
926  void output(std::ostream &outfile);
927 
928  /// Output at f(n_plot) points in each element
929  void output(std::ostream &outfile, const unsigned &n_plot);
930 
931  /// Output for all elements (C-style output)
932  void output(FILE* file_pt);
933 
934  /// Output at f(n_plot) points in each element (C-style output)
935  void output(FILE* file_pt, const unsigned &nplot);
936 
937  /// Output for all elements
938  void output(const std::string& output_filename)
939  {
940  std::ofstream outfile;
941  outfile.open(output_filename.c_str());
942  output(outfile);
943  outfile.close();
944  }
945 
946  /// Output at f(n_plot) points in each element
947  void output(const std::string& output_filename, const unsigned &n_plot)
948  {
949  std::ofstream outfile;
950  outfile.open(output_filename.c_str());
951  output(outfile,n_plot);
952  outfile.close();
953  }
954 
955  /// Output a given Vector function at f(n_plot) points in each element
956  void output_fct(std::ostream &outfile, const unsigned &n_plot,
958 
959  /// \short Output a given time-dep. Vector function at f(n_plot) points in
960  /// each element
961  void output_fct(std::ostream &outfile, const unsigned &n_plot,
962  const double& time,
964 
965  /// Output the nodes on the boundaries (into separate tecplot zones)
966  void output_boundaries(std::ostream &outfile);
967 
968  /// Output the nodes on the boundaries (into separate tecplot zones).
969  /// Specify filename
970  void output_boundaries(const std::string& output_filename)
971  {
972  std::ofstream outfile;
973  outfile.open(output_filename.c_str());
974  output_boundaries(outfile);
975  outfile.close();
976  }
977 
978  /// Output the nodes on the boundary and their respective boundary
979  /// coordinates(into separate tecplot zones)
980  virtual void output_boundary_coordinates(const unsigned &b,
981  std::ostream &outfile)
982  {
983  std::ostringstream error_stream;
984  error_stream
985  << "Empty default output_boundary_coordinates() method called.\n"
986  << "This should be overloaded on the specific mesh that you are using\n";
987  throw OomphLibError(error_stream.str(),
988  "Mesh::output_boundary_coordinates",
989  OOMPH_EXCEPTION_LOCATION);
990  }
991 
992  /// Output the nodes on the boundary and their respective boundary
993  /// coordinates(into separate tecplot zones)
994  void output_boundary_coordinates(const unsigned &b,
995  const std::string& output_filename)
996  {
997  std::ofstream outfile;
998  outfile.open(output_filename.c_str());
999  output_boundary_coordinates(b, outfile);
1000  outfile.close();
1001  }
1002 
1003  /// Output the nodes on the boundaries and their respective boundary
1004  /// coordinates(into separate tecplot zones)
1005  void output_boundaries_coordinates(std::ostream &outfile);
1006 
1007  /// Output the nodes on the boundaries and their respective boundary
1008  /// coordinates(into separate tecplot zones)
1009  void output_boundaries_coordinates(const std::string& output_filename)
1010  {
1011  std::ofstream outfile;
1012  outfile.open(output_filename.c_str());
1014  outfile.close();
1015  }
1016 
1017  /// \short Assign initial values for an impulsive start
1019 
1020  /// \short Shift time-dependent data along for next timestep:
1021  /// Deal with nodal Data/positions and the element's internal
1022  /// Data
1023  void shift_time_values();
1024 
1025 
1026  /// \short Calculate predictions for all Data and positions associated
1027  /// with the mesh, usually used in adaptive time-stepping.
1028  void calculate_predictions();
1029 
1030  /// \short Set the timestepper associated with all nodal and elemental
1031  /// data stored in the mesh.
1033  const bool &preserve_existing_data)
1034  {
1035  this->set_nodal_time_stepper(time_stepper_pt,preserve_existing_data);
1036  this->set_elemental_internal_time_stepper(time_stepper_pt,
1037  preserve_existing_data);
1038  }
1039 
1040  /// \short Function that can be used to set any additional timestepper data
1041  /// stored at the Mesh (as opposed to nodal and elemental) levels. This
1042  /// is virtual so that it can be overloaded in the appropriate Meshes.
1043  /// Examples include the SpineMeshes and adaptive triangle and tet meshes
1044  virtual void set_mesh_level_time_stepper(TimeStepper* const &time_stepper_pt,
1045  const bool &preserve_existing_data);
1046 
1047  /// \short Set consistent values for pinned data in continuation
1049  ContinuationStorageScheme* const &continuation_stepper_pt);
1050 
1051  /// \short Does the double pointer correspond to any mesh data
1052  bool does_pointer_correspond_to_mesh_data(double* const &parameter_pt);
1053 
1054  /// \short Set the timestepper associated with the nodal data in the mesh
1055  void set_nodal_time_stepper(TimeStepper* const &time_stepper_pt,
1056  const bool &preserve_existing_data);
1057 
1058  /// \short Set the timestepper associated with the internal data stored
1059  /// within elements in the meah
1060  void set_elemental_internal_time_stepper(TimeStepper* const &time_stepper_pt,
1061  const bool &preserve_existing_data);
1062 
1063  /// \short Compute norm of solution by summing contributions of
1064  /// compute_norm(...) for all constituent elements in the mesh.
1065  /// What that norm means depends on what's defined in the element's
1066  /// function; may need to take the square root afterwards if the elements
1067  /// compute the square of the L2 norm, say.
1068  virtual void compute_norm(double& norm)
1069  {
1070  //Initialse the norm
1071  norm=0.0;
1072 
1073  //Per-element norm
1074  double el_norm=0;
1075 
1076  //Loop over the elements
1077  unsigned long n_element = Element_pt.size();
1078  for(unsigned long e=0;e<n_element;e++)
1079  {
1081 
1082 #ifdef OOMPH_HAS_MPI
1083  //Compute error for each non-halo element
1084  if (!(el_pt->is_halo()))
1085 #endif
1086  {
1087  el_pt->compute_norm(el_norm);
1088  }
1089  norm+=el_norm;
1090  }
1091  }
1092 
1093 
1094  /// \short Plot error when compared against a given exact solution.
1095  /// Also returns the norm of the error and that of the exact solution
1096  virtual void compute_error(std::ostream &outfile,
1098  const double& time,
1099  double& error, double& norm)
1100  {
1101  //Initialse the norm and error
1102  norm=0.0; error=0.0;
1103  //Per-element norm and error
1104  double el_error,el_norm;
1105 
1106  //Loop over the elements
1107  unsigned long Element_pt_range = Element_pt.size();
1108  for(unsigned long e=0;e<Element_pt_range;e++)
1109  {
1110  // Try to cast to FiniteElement
1111  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
1112  if (el_pt==0)
1113  {
1114  throw OomphLibError(
1115  "Can't execute compute_error(...) for non FiniteElements",
1116  OOMPH_CURRENT_FUNCTION,
1117  OOMPH_EXCEPTION_LOCATION);
1118  }
1119 
1120  // Reset elemental errors and norms
1121  el_norm=0.0;
1122  el_error=0.0;
1123 #ifdef OOMPH_HAS_MPI
1124  //Compute error for each non-halo element
1125  if (!(el_pt->is_halo()))
1126 #endif
1127  {
1128  el_pt->compute_error(outfile,exact_soln_pt,time,el_error,el_norm);
1129  }
1130  //Add each element error to the global errors
1131  norm+=el_norm; error+=el_error;
1132  }
1133  }
1134 
1135  /// \short Plot error when compared against a given time-depdendent
1136  /// exact solution. Also returns the norm of the error and
1137  /// that of the exact solution
1138  virtual void compute_error(std::ostream &outfile,
1140  double& error, double& norm)
1141  {
1142  //Initialise norm and error
1143  norm=0.0; error=0.0;
1144  //Per-element norm and error
1145  double el_error,el_norm;
1146 
1147  //Loop over the elements
1148  unsigned long Element_pt_range = Element_pt.size();
1149  for(unsigned long e=0;e<Element_pt_range;e++)
1150  {
1151  // Try to cast to FiniteElement
1152  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
1153  if (el_pt==0)
1154  {
1155  throw OomphLibError(
1156  "Can't execute compute_error(...) for non FiniteElements",
1157  OOMPH_CURRENT_FUNCTION,
1158  OOMPH_EXCEPTION_LOCATION);
1159  }
1160  // Reset elemental errors and norms
1161  el_norm=0.0;
1162  el_error=0.0;
1163  //Calculate the elemental errors for each non-halo element
1164 #ifdef OOMPH_HAS_MPI
1165  if (!(el_pt->is_halo()))
1166 #endif
1167  {
1168  el_pt->compute_error(outfile,exact_soln_pt,el_error,el_norm);
1169  }
1170  //Add each elemental error to the global error
1171  norm+=el_norm; error+=el_error;
1172  }
1173  }
1174 
1175  /// \short Plot error when compared against a given time-depdendent
1176  /// exact solution. Also returns the norm of the error and
1177  /// that of the exact solution. Version with vectors of norms and errors so
1178  /// that different variables' norms and errors can be returned individually
1179  virtual void compute_error(
1180  std::ostream &outfile,
1182  const double& time,
1183  Vector<double>& error, Vector<double>& norm)
1184  {
1185  //Initialise norm and error
1186  unsigned n_error=error.size();
1187  unsigned n_norm=norm.size();
1188  for(unsigned i=0;i<n_error;i++)
1189  {
1190  error[i]=0.0;
1191  }
1192  for(unsigned i=0;i<n_norm;i++)
1193  {
1194  norm[i]=0.0;
1195  }
1196  //Per-element norm and error
1197  Vector<double> el_error(n_error),el_norm(n_norm);
1198 
1199  //Loop over the elements
1200  unsigned long Element_pt_range = Element_pt.size();
1201  for(unsigned long e=0;e<Element_pt_range;e++)
1202  {
1203  // Try to cast to FiniteElement
1204  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
1205  if (el_pt==0)
1206  {
1207  throw OomphLibError(
1208  "Can't execute compute_error(...) for non FiniteElements",
1209  OOMPH_CURRENT_FUNCTION,
1210  OOMPH_EXCEPTION_LOCATION);
1211  }
1212  // Reset elemental errors and norms
1213  for(unsigned i=0;i<n_error;i++)
1214  {
1215  el_error[i]=0.0;
1216  }
1217  for(unsigned i=0;i<n_norm;i++)
1218  {
1219  el_norm[i]=0.0;
1220  }
1221  //Calculate the elemental errors for each non-halo element
1222 #ifdef OOMPH_HAS_MPI
1223  if (!(el_pt->is_halo()))
1224 #endif
1225  {
1226  el_pt->compute_error(outfile,exact_soln_pt,time,el_error,el_norm);
1227  }
1228  //Add each elemental error to the global error
1229  for(unsigned i=0;i<n_error;i++)
1230  {
1231  error[i]+=el_error[i];
1232  }
1233  for(unsigned i=0;i<n_norm;i++)
1234  {
1235  norm[i]+=el_norm[i];
1236  }
1237  }
1238  }
1239 
1240  /// \short Plot error when compared against a given time-depdendent
1241  /// exact solution. Also returns the norm of the error and
1242  /// that of the exact solution. Version with vectors of norms and errors so
1243  /// that different variables' norms and errors can be returned individually
1244  virtual void compute_error(std::ostream &outfile,
1246  Vector<double>& error, Vector<double>& norm)
1247  {
1248  //Initialise norm and error
1249  unsigned n_error=error.size();
1250  unsigned n_norm=norm.size();
1251  for(unsigned i=0;i<n_error;i++)
1252  {
1253  error[i]=0.0;
1254  }
1255  for(unsigned i=0;i<n_norm;i++)
1256  {
1257  norm[i]=0.0;
1258  }
1259  //Per-element norm and error
1260  Vector<double> el_error(n_error),el_norm(n_norm);
1261 
1262  //Loop over the elements
1263  unsigned long Element_pt_range = Element_pt.size();
1264  for(unsigned long e=0;e<Element_pt_range;e++)
1265  {
1266  // Try to cast to FiniteElement
1267  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
1268  if (el_pt==0)
1269  {
1270  throw OomphLibError(
1271  "Can't execute compute_error(...) for non FiniteElements",
1272  OOMPH_CURRENT_FUNCTION,
1273  OOMPH_EXCEPTION_LOCATION);
1274  }
1275  // Reset elemental errors and norms
1276  for(unsigned i=0;i<n_error;i++)
1277  {
1278  el_error[i]=0.0;
1279  }
1280  for(unsigned i=0;i<n_norm;i++)
1281  {
1282  el_norm[i]=0.0;
1283  }
1284  //Calculate the elemental errors for each non-halo element
1285 #ifdef OOMPH_HAS_MPI
1286  if (!(el_pt->is_halo()))
1287 #endif
1288  {
1289  el_pt->compute_error(outfile,exact_soln_pt,el_error,el_norm);
1290  }
1291  //Add each elemental error to the global error
1292  for(unsigned i=0;i<n_error;i++)
1293  {
1294  error[i]+=el_error[i];
1295  }
1296  for(unsigned i=0;i<n_norm;i++)
1297  {
1298  norm[i]+=el_norm[i];
1299  }
1300  }
1301  }
1302 
1303 
1304  /// Boolean to indicate if Mesh has been distributed
1305  bool is_mesh_distributed() const
1306  {
1307 #ifdef OOMPH_HAS_MPI
1308  return communicator_pt() != 0;
1309 #else
1310  // Can't be distributed without mpi
1311  return false;
1312 #endif
1313  }
1314 
1315  /// Read-only access fct to communicator (Null if mesh is not distributed,
1316  /// i.e. if we don't have mpi).
1318  {
1319 #ifdef OOMPH_HAS_MPI
1320  return Comm_pt;
1321 #else
1322  return 0;
1323 #endif
1324  }
1325 
1326 
1327 #ifdef OOMPH_HAS_MPI
1328 
1329  /// Function to set communicator (mesh is assumed to be distributed if the
1330  /// communicator pointer is non-null). Only defined if mpi is enabled
1331  /// becaus Comm_pt itself is only defined when mpi is enabled.
1332  void set_communicator_pt(OomphCommunicator* comm_pt) {Comm_pt = comm_pt;}
1333 
1334  /// \short Call this function to keep all the elements as halo elements
1336 
1337  /// \short Calll this function to unset the flag that keeps all elements
1338  /// in the mesh as halo elements
1340 
1341  /// \short Distribute the problem and doc; make this virtual to allow
1342  /// overloading for particular meshes where further work is required.
1343  /// Add to vector of pointers to deleted elements.
1344  virtual void distribute(OomphCommunicator* comm_pt,
1345  const Vector<unsigned>& element_domain,
1346  Vector<GeneralisedElement*>& deleted_element_pt,
1347  DocInfo& doc_info,
1348  const bool& report_stats,
1349  const bool& overrule_keep_as_halo_element_status);
1350 
1351  /// \short Distribute the problem
1352  /// Add to vector of pointers to deleted elements.
1354  const Vector<unsigned>& element_domain,
1355  Vector<GeneralisedElement*>& deleted_element_pt,
1356  const bool& report_stats=false)
1357  {
1358  DocInfo doc_info;
1359  doc_info.disable_doc();
1360  bool overrule_keep_as_halo_element_status=false;
1361  return distribute(comm_pt,element_domain,deleted_element_pt,
1362  doc_info,report_stats,
1363  overrule_keep_as_halo_element_status);
1364  }
1365 
1366  /// \short (Irreversibly) prune halo(ed) elements and nodes, usually
1367  /// after another round of refinement, to get rid of
1368  /// excessively wide halo layers. Note that the current
1369  /// mesh will be now regarded as the base mesh and no unrefinement
1370  /// relative to it will be possible once this function
1371  /// has been called.
1373  deleted_element_pt,
1374  const bool& report_stats=false)
1375  {
1376  DocInfo doc_info;
1377  doc_info.disable_doc();
1378  prune_halo_elements_and_nodes(deleted_element_pt,
1379  doc_info,report_stats);
1380  }
1381 
1382 
1383  /// \short (Irreversibly) prune halo(ed) elements and nodes, usually
1384  /// after another round of refinement, to get rid of
1385  /// excessively wide halo layers. Note that the current
1386  /// mesh will be now regarded as the base mesh and no unrefinement
1387  /// relative to it will be possible once this function
1388  /// has been called.
1390  deleted_element_pt,
1391  DocInfo& doc_info,
1392  const bool& report_stats);
1393 
1394  /// \short Get efficiency of mesh distribution: In an ideal distribution
1395  /// without halo overhead, each processor would only hold its own
1396  /// elements. Efficieny per processor = (number of non-halo elements)/
1397  /// (total number of elements).
1398  void get_efficiency_of_mesh_distribution(double& av_efficiency,
1399  double& max_efficiency,
1400  double& min_efficiency);
1401 
1402  /// Doc the mesh distribution, to be processed with tecplot macros
1403  void doc_mesh_distribution(DocInfo& doc_info);
1404 
1405  /// Check halo and shared schemes on the mesh
1406  void check_halo_schemes(DocInfo& doc_info,
1407  double& max_permitted_error_for_halo_check);
1408 
1409  /// \short Classify the halo and haloed nodes in the mesh. Virtual
1410  /// so it can be overloaded to perform additional functionality
1411  /// (such as synchronising hanging nodes) in refineable meshes, say.
1412  virtual void classify_halo_and_haloed_nodes(DocInfo& doc_info,
1413  const bool& report_stats);
1414 
1415  /// Classify the halo and haloed nodes in the mesh. Virtual
1416  /// so it can be overloaded to perform additional functionality
1417  /// (such as synchronising hanging nodes) in refineable meshes, say.
1418  virtual void classify_halo_and_haloed_nodes(const bool& report_stats=false)
1419  {
1420  DocInfo doc_info;
1421  doc_info.disable_doc();
1422  classify_halo_and_haloed_nodes(doc_info,report_stats);
1423  }
1424 
1425  /// \short Synchronise shared node lookup schemes to cater for the
1426  /// the case where:
1427  /// (1) a certain node on the current processor is halo with proc p
1428  /// (i.e. its non-halo counterpart lives on processor p)
1429  /// (2) that node is also exists (also as a halo) on another processor
1430  /// (q, say) where its non-halo counter part is also known to be
1431  /// on processor p.
1432  /// However, without calling this function the current processor does not
1433  /// necessarily know that it shares a node with processor q. This
1434  /// information can be required, e.g. when synchronising hanging node
1435  /// schemes over all processors.
1436  void synchronise_shared_nodes(const bool& report_stats);
1437 
1438 
1439  /// \short Get all the halo data stored in the mesh and add pointers to
1440  /// the data to the map, indexed by global equation number
1441  void get_all_halo_data(std::map<unsigned,double*> &map_of_halo_data);
1442 
1443  /// \short Return vector of halo elements in this Mesh
1444  /// whose non-halo counterpart is held on processor p.
1446  {
1447  // Prepare vector
1448  Vector<GeneralisedElement*> vec_el_pt;
1449 
1450  // Loop over all root halo elements
1451  unsigned nelem=nroot_halo_element(p);
1452  for (unsigned e=0;e<nelem;e++)
1453  {
1455 
1456  // Is it a refineable element?
1457  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(el_pt);
1458  if (ref_el_pt!=0)
1459  {
1460  // Vector of pointers to leaves in tree emanating from
1461  // current root halo element
1462  Vector<Tree*> leaf_pt;
1463  ref_el_pt->tree_pt()->stick_leaves_into_vector(leaf_pt);
1464 
1465  // Loop over leaves and add their objects (the finite elements)
1466  // to vector
1467  unsigned nleaf=leaf_pt.size();
1468  for (unsigned l=0;l<nleaf;l++)
1469  {
1470  vec_el_pt.push_back(leaf_pt[l]->object_pt());
1471  }
1472  }
1473  else
1474  {
1475  vec_el_pt.push_back(el_pt);
1476  }
1477  }
1478  return vec_el_pt;
1479  }
1480 
1481 
1482  /// \short Return vector of haloed elements in this Mesh
1483  /// whose haloing counterpart is held on processor p.
1485  {
1486  // Prepare vector
1487  Vector<GeneralisedElement*> vec_el_pt;
1488 
1489  // Loop over all root haloed elements
1490  unsigned nelem=nroot_haloed_element(p);
1491  for (unsigned e=0;e<nelem;e++)
1492  {
1494 
1495  // Is it a refineable element?
1496  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(el_pt);
1497  if (ref_el_pt!=0)
1498  {
1499  // Vector of pointers to leaves in tree emanating from
1500  // current root haloed element
1501  Vector<Tree*> leaf_pt;
1502  ref_el_pt->tree_pt()->stick_leaves_into_vector(leaf_pt);
1503 
1504  // Loop over leaves and add their objects (the finite elements)
1505  // to vector
1506  unsigned nleaf=leaf_pt.size();
1507  for (unsigned l=0;l<nleaf;l++)
1508  {
1509  vec_el_pt.push_back(leaf_pt[l]->object_pt());
1510  }
1511  }
1512  else
1513  {
1514  vec_el_pt.push_back(el_pt);
1515  }
1516  }
1517  return vec_el_pt;
1518  }
1519 
1520  /// \short Total number of non-halo elements in this mesh (Costly call computes
1521  /// result on the fly)
1523  {
1524  unsigned count=0;
1525  unsigned n=nelement();
1526  for (unsigned e=0;e<n;e++)
1527  {
1528  if (!(element_pt(e)->is_halo())) count++;
1529  }
1530  return count;
1531  }
1532 
1533 
1534  /// \short Total number of root halo elements in this Mesh
1536  {
1537  unsigned n=0;
1538  for (std::map<unsigned,Vector<GeneralisedElement*> >::iterator it=
1539  Root_halo_element_pt.begin();it!=Root_halo_element_pt.end();it++)
1540  {
1541  n+=it->second.size();
1542  }
1543  return n;
1544  }
1545 
1546 
1547  /// \short Number of root halo elements in this Mesh whose non-halo
1548  /// counterpart is held on processor p.
1549  unsigned nroot_halo_element(const unsigned& p)
1550  {
1551  return Root_halo_element_pt[p].size();
1552  }
1553 
1554 
1555  /// \short Vector of pointers to root halo elements in this Mesh
1556  /// whose non-halo counterpart is held on processor p.
1558  {
1559  return Root_halo_element_pt[p];
1560  }
1561 
1562 
1563  /// \short Access fct to the e-th root halo element in this Mesh
1564  /// whose non-halo counterpart is held on processor p.
1566  const unsigned& e)
1567  {
1568  return Root_halo_element_pt[p][e];
1569  }
1570 
1571 
1572  /// \short Add root halo element whose non-halo counterpart is held
1573  /// on processor p to this Mesh.
1574  void add_root_halo_element_pt(const unsigned& p, GeneralisedElement*& el_pt)
1575  {
1576  Root_halo_element_pt[p].push_back(el_pt);
1577  el_pt->set_halo(p);
1578  }
1579 
1580  /// \short Total number of halo nodes in this Mesh
1581  unsigned nhalo_node()
1582  {
1583  unsigned n=0;
1584  for (std::map<unsigned,Vector<Node*> >::iterator it=
1585  Halo_node_pt.begin();it!=Halo_node_pt.end();it++)
1586  {
1587  n+=it->second.size();
1588  }
1589  return n;
1590  }
1591 
1592  /// \short Number of halo nodes in this Mesh whose non-halo counterpart
1593  /// is held on processor p.
1594  unsigned nhalo_node(const unsigned& p)
1595  {
1596  //Memory saving version of: return Halo_node_pt[p].size();
1597  std::map<unsigned, Vector<Node*> >::iterator it=
1598  Halo_node_pt.find(p);
1599  if (it==Halo_node_pt.end())
1600  {
1601  return 0;
1602  }
1603  return (*it).second.size();
1604  }
1605 
1606  /// \short Add halo node whose non-halo counterpart is held
1607  /// on processor p to the storage scheme for halo nodes.
1608  void add_halo_node_pt(const unsigned& p, Node*& nod_pt)
1609  {
1610  Halo_node_pt[p].push_back(nod_pt);
1611  }
1612 
1613 
1614  /// \short Access fct to the j-th halo node in this Mesh
1615  /// whose non-halo counterpart is held on processor p.
1616  Node* halo_node_pt(const unsigned& p, const unsigned& j)
1617  {
1618  return Halo_node_pt[p][j];
1619  }
1620 
1621 
1622  /// \short Total number of root haloed elements in this Mesh
1624  {
1625  unsigned n=0;
1626  for (std::map<unsigned,Vector<GeneralisedElement*> >::iterator it=
1627  Root_haloed_element_pt.begin();it!=Root_haloed_element_pt.end();it++)
1628  {
1629  n+=it->second.size();
1630  }
1631  return n;
1632  }
1633 
1634  /// \short Number of root haloed elements in this Mesh whose non-halo
1635  /// counterpart is held on processor p.
1636  unsigned nroot_haloed_element(const unsigned& p)
1637  {
1638  //Memory saving version of: return Root_haloed_element_pt[p].size();
1639  std::map<unsigned, Vector<GeneralisedElement*> >::iterator it=
1640  Root_haloed_element_pt.find(p);
1641  if (it==Root_haloed_element_pt.end())
1642  {
1643  return 0;
1644  }
1645  return (*it).second.size();
1646  }
1647 
1648 
1649  /// \short Vector of pointers to root haloed elements in this Mesh
1650  /// whose non-halo counterpart is held on processor p.
1652  {
1653  //Memory saving version of: return Root_haloed_element_pt[p];
1654  std::map<unsigned, Vector<GeneralisedElement*> >::iterator it=
1655  Root_haloed_element_pt.find(p);
1656  if (it==Root_haloed_element_pt.end())
1657  {
1659  return tmp;
1660  }
1661  return (*it).second;
1662  }
1663 
1664  /// \short Access fct to the e-th root haloed element in this Mesh
1665  /// whose non-halo counterpart is held on processor p.
1667  const unsigned& e)
1668  {
1669  return Root_haloed_element_pt[p][e];
1670  }
1671 
1672  /// \short Add root haloed element whose non-halo counterpart is held
1673  /// on processor p to the storage scheme for haloed elements.
1674  /// Note: This does not add the element to the storage scheme
1675  /// for elements as it's understood to naturally live on this
1676  /// processor anyway!
1677  void add_root_haloed_element_pt(const unsigned& p, GeneralisedElement*& el_pt)
1678  {
1679  Root_haloed_element_pt[p].push_back(el_pt);
1680  }
1681 
1682 
1683  /// \short Total number of haloed nodes in this Mesh
1684  unsigned nhaloed_node()
1685  {
1686  unsigned n=0;
1687  for (std::map<unsigned,Vector<Node*> >::iterator it=
1688  Haloed_node_pt.begin();it!=Haloed_node_pt.end();it++)
1689  {
1690  n+=it->second.size();
1691  }
1692  return n;
1693  }
1694 
1695 
1696  /// \short Number of haloed nodes in this Mesh whose haloed counterpart
1697  /// is held on processor p.
1698  unsigned nhaloed_node(const unsigned& p)
1699  {
1700  // Memory saving version of: return Haloed_node_pt[p].size();
1701  std::map<unsigned, Vector<Node*> >::iterator it=
1702  Haloed_node_pt.find(p);
1703  if (it==Haloed_node_pt.end())
1704  {
1705  return 0;
1706  }
1707  return (*it).second.size();
1708  }
1709 
1710  /// \short Access fct to the j-th haloed node in this Mesh
1711  /// whose halo counterpart is held on processor p.
1712  Node* haloed_node_pt(const unsigned& p, const unsigned& j)
1713  {
1714  return Haloed_node_pt[p][j];
1715  }
1716 
1717  /// \short Add haloed node whose halo counterpart is held
1718  /// on processor p to the storage scheme for haloed nodes.
1719  void add_haloed_node_pt(const unsigned& p, Node*& nod_pt)
1720  {
1721  Haloed_node_pt[p].push_back(nod_pt);
1722  }
1723 
1724  /// Bool for output of halo elements
1726 
1727 
1728  /// \short Function to suppress resizing of halo nodes -- optmisation
1729  /// but call it at your own risk!
1731  {
1733  }
1734 
1735 
1736  /// \short Function to (re-)enable resizing of halo nodes -- this returns
1737  /// things to the default behaviour.
1739  {
1741  }
1742 
1743  /// Function to disable halo element output
1745  {
1746  Output_halo_elements=false;
1747  }
1748 
1749  /// Function to enable halo element output
1751  {
1752  Output_halo_elements=true;
1753  }
1754 
1755  /// \short Total number of shared nodes in this Mesh
1756  unsigned nshared_node()
1757  {
1758  unsigned n=0;
1759  for (std::map<unsigned,Vector<Node*> >::iterator it=
1760  Shared_node_pt.begin();it!=Shared_node_pt.end();it++)
1761  {
1762  n+=it->second.size();
1763  }
1764  return n;
1765  }
1766 
1767  /// \short Doc shared nodes
1769  {
1770  for (std::map<unsigned,Vector<Node*> >::iterator it=
1771  Shared_node_pt.begin();it!=Shared_node_pt.end();it++)
1772  {
1773  unsigned n=it->second.size();
1774  for (unsigned j=0;j<n;j++)
1775  {
1776  oomph_info << "Shared node with proc " << it->first << " ";
1777  Node* nod_pt=it->second[j];
1778  unsigned ndim=nod_pt->ndim();
1779  for (unsigned i=0;i<ndim;i++)
1780  {
1781  oomph_info << nod_pt->x(i) << " ";
1782  }
1783  oomph_info << nod_pt->is_hanging() << " ";
1784  oomph_info << j << " " << nod_pt << std::endl;
1785  }
1786  }
1787  }
1788 
1789  /// \short Number of shared nodes in this Mesh who have a counterpart
1790  /// on processor p.
1791  unsigned nshared_node(const unsigned& p)
1792  {
1793  // Memory saving version of: return Shared_node_pt[p].size();
1794  std::map<unsigned, Vector<Node*> >::iterator it=
1795  Shared_node_pt.find(p);
1796  if (it==Shared_node_pt.end())
1797  {
1798  return 0;
1799  }
1800  return (*it).second.size();
1801  }
1802 
1803  /// \short Access fct to the j-th shared node in this Mesh
1804  /// who has a counterpart on processor p.
1805  Node* shared_node_pt(const unsigned& p, const unsigned& j)
1806  {
1807  return Shared_node_pt[p][j];
1808  }
1809 
1810  /// \short Get vector of pointers to shared nodes with processor p.
1811  /// Required for faster search in
1812  /// Missing_masters_functions::add_external_haloed_node_helper() and
1813  /// Missing_masters_functions::add_external_haloed_master_node_helper()
1815  {
1816  unsigned np=nshared_node(p);
1817  shared_node_pt.resize(np);
1818  for (unsigned j=0;j<np;j++)
1819  {
1820  shared_node_pt[j] = Shared_node_pt[p][j];
1821  }
1822  }
1823 
1824 
1825  /// \short Add shared node whose counterpart is held
1826  /// on processor p to the storage scheme for shared nodes.
1827  /// (NB: ensure that this routine is called twice, once for each process)
1828  void add_shared_node_pt(const unsigned& p, Node*& nod_pt)
1829  {
1830  Shared_node_pt[p].push_back(nod_pt);
1831  }
1832 
1833  /// \short Get halo node stats for this distributed mesh:
1834  /// Average/max/min number of halo nodes over all processors.
1835  /// \b Careful: Involves MPI Broadcasts and must therefore
1836  /// be called on all processors!
1837  void get_halo_node_stats(double& av_number,
1838  unsigned& max_number,
1839  unsigned& min_number);
1840 
1841  /// \short Get haloed node stats for this distributed mesh:
1842  /// Average/max/min number of haloed nodes over all processors.
1843  /// \b Careful: Involves MPI Broadcasts and must therefore
1844  /// be called on all processors!
1845  void get_haloed_node_stats(double& av_number,
1846  unsigned& max_number,
1847  unsigned& min_number);
1848 
1849  // External halo(ed) elements are "source/other" elements which are
1850  // on different processes to the element for which they are the source
1851 
1852  /// \short Output all external halo elements
1853  void output_external_halo_elements(std::ostream &outfile,
1854  const unsigned &n_plot=5)
1855  {
1856  for (std::map<unsigned,Vector<GeneralisedElement*> >::iterator it=
1857  External_halo_element_pt.begin();
1858  it!=External_halo_element_pt.end();it++)
1859  {
1860  unsigned p=(*it).first;
1861  output_external_halo_elements(p,outfile,n_plot);
1862  }
1863  }
1864 
1865  /// \short Output all external halo elements with processor p
1866  void output_external_halo_elements(const unsigned& p, std::ostream &outfile,
1867  const unsigned &n_plot=5)
1868  {
1869  unsigned nel=External_halo_element_pt[p].size();
1870  for (unsigned e=0;e<nel;e++)
1871  {
1872  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(
1874  if (fe_pt!=0)
1875  {
1876  fe_pt->output(outfile,n_plot);
1877  }
1878  }
1879  }
1880 
1881 
1882  /// \short Output all external haloed elements
1883  void output_external_haloed_elements(std::ostream &outfile,
1884  const unsigned &n_plot=5)
1885  {
1886  for (std::map<unsigned,Vector<GeneralisedElement*> >::iterator it=
1888  it!=External_haloed_element_pt.end();it++)
1889  {
1890  unsigned p=(*it).first;
1891  output_external_haloed_elements(p,outfile,n_plot);
1892  }
1893  }
1894 
1895  /// \short Output all external haloed elements with processor p
1896  void output_external_haloed_elements(const unsigned& p, std::ostream &outfile,
1897  const unsigned &n_plot=5)
1898  {
1899  unsigned nel=External_haloed_element_pt[p].size();
1900  for (unsigned e=0;e<nel;e++)
1901  {
1902  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(
1904  if (fe_pt!=0)
1905  {
1906  fe_pt->output(outfile,n_plot);
1907  }
1908  }
1909  }
1910 
1911 
1912  /// \short Total number of external halo elements in this Mesh
1914  {
1915  unsigned n=0;
1916  for (std::map<unsigned,Vector<GeneralisedElement*> >::iterator it=
1917  External_halo_element_pt.begin();
1918  it!=External_halo_element_pt.end();it++)
1919  {
1920  n+=it->second.size();
1921  }
1922  return n;
1923  }
1924 
1925  /// \short Number of external halo elements in this Mesh whose non-halo
1926  /// counterpart is held on processor p.
1927  unsigned nexternal_halo_element(const unsigned& p)
1928  {
1929  // Memory saving version of: return External_halo_element_pt[p].size();
1930  std::map<unsigned, Vector<GeneralisedElement*> >::iterator it=
1931  External_halo_element_pt.find(p);
1932  if (it==External_halo_element_pt.end())
1933  {
1934  return 0;
1935  }
1936  return (*it).second.size();
1937  }
1938 
1939  /// \short Access fct to the e-th external halo element in this Mesh
1940  /// whose non-halo counterpart is held on processor p.
1942  const unsigned& e)
1943  {
1944  return External_halo_element_pt[p][e];
1945  }
1946 
1947  /// \short Add external halo element whose non-halo counterpart is held
1948  /// on processor p to this Mesh.
1949  void add_external_halo_element_pt(const unsigned& p,
1950  GeneralisedElement*& el_pt)
1951  {
1952  External_halo_element_pt[p].push_back(el_pt);
1953  el_pt->set_halo(p);
1954  }
1955 
1956  /// \short Total number of external haloed elements in this Mesh
1958  {
1959  unsigned n=0;
1960  for (std::map<unsigned,Vector<GeneralisedElement*> >::iterator it=
1962  it!=External_haloed_element_pt.end();it++)
1963  {
1964  n+=it->second.size();
1965  }
1966  return n;
1967  }
1968 
1969  /// \short Number of external haloed elements in this Mesh whose non-halo
1970  /// counterpart is held on processor p.
1971  unsigned nexternal_haloed_element(const unsigned& p)
1972  {
1973  // Memory saving version of: return External_haloed_element_pt[p].size();
1974  std::map<unsigned, Vector<GeneralisedElement*> >::iterator it=
1976  if (it==External_haloed_element_pt.end())
1977  {
1978  return 0;
1979  }
1980  return (*it).second.size();
1981  }
1982 
1983  /// \short Access fct to the e-th external haloed element in this Mesh
1984  /// whose non-halo counterpart is held on processor p.
1986  const unsigned& e)
1987  {
1988  return External_haloed_element_pt[p][e];
1989  }
1990 
1991  /// \short Add external haloed element whose non-halo counterpart is held
1992  /// on processor p to the storage scheme for haloed elements.
1993  unsigned add_external_haloed_element_pt(const unsigned& p,
1994  GeneralisedElement*& el_pt);
1995 
1996  /// \short Total number of external halo nodes in this Mesh
1998  {
1999  unsigned n=0;
2000  for (std::map<unsigned,Vector<Node*> >::iterator it=
2001  External_halo_node_pt.begin();it!=External_halo_node_pt.end();it++)
2002  {
2003  n+=it->second.size();
2004  }
2005  return n;
2006  }
2007 
2008 
2009  /// \short Get vector of pointers to all external halo nodes
2011  {
2012  unsigned n_total=nexternal_halo_node();
2013  external_halo_node_pt.reserve(n_total);
2014  for (std::map<unsigned,Vector<Node*> >::iterator it=
2015  External_halo_node_pt.begin();it!=External_halo_node_pt.end();it++)
2016  {
2017  unsigned np=(it->second).size();
2018  for (unsigned j=0;j<np;j++)
2019  {
2020  external_halo_node_pt.push_back((it->second)[j]);
2021  }
2022  }
2023 #ifdef PARANOID
2024  if (external_halo_node_pt.size()!=n_total)
2025  {
2026  std::ostringstream error_stream;
2027  error_stream
2028  << "Total number of external halo nodes, "
2029  << n_total << " doesn't match number of entries \n in vector, "
2030  << external_halo_node_pt.size() << std::endl;
2031  throw OomphLibError(error_stream.str(),
2032  OOMPH_CURRENT_FUNCTION,
2033  OOMPH_EXCEPTION_LOCATION);
2034  }
2035 #endif
2036  }
2037 
2038 
2039  /// \short Number of external halo nodes in this Mesh whose non-halo
2040  /// (external) counterpart is held on processor p.
2041  unsigned nexternal_halo_node(const unsigned& p)
2042  {
2043  // Memory saving version of: return External_halo_node_pt[p].size();
2044  std::map<unsigned, Vector<Node*> >::iterator it=
2045  External_halo_node_pt.find(p);
2046  if (it==External_halo_node_pt.end())
2047  {
2048  return 0;
2049  }
2050  return (*it).second.size();
2051  }
2052 
2053  /// \short Add external halo node whose non-halo (external) counterpart
2054  /// is held on processor p to the storage scheme for halo nodes.
2055  void add_external_halo_node_pt(const unsigned& p, Node*& nod_pt)
2056  {
2057  External_halo_node_pt[p].push_back(nod_pt);
2058  nod_pt->set_halo(p);
2059  }
2060 
2061 
2062  /// \short Access fct to the j-th external halo node in this Mesh
2063  /// whose non-halo external counterpart is held on processor p.
2064  Node* &external_halo_node_pt(const unsigned& p, const unsigned& j)
2065  {
2066  return External_halo_node_pt[p][j];
2067  }
2068 
2069  /// \short Access fct to vector of external halo node in this Mesh
2070  /// whose non-halo external counterpart is held on processor p. (read only)
2072  {
2073  std::map<unsigned, Vector<Node*> >::iterator it=
2074  External_halo_node_pt.find(p);
2075  if (it==External_halo_node_pt.end())
2076  {
2077  Vector<Node*> tmp;
2078  return tmp;
2079  }
2080  return (*it).second;
2081  }
2082 
2083  /// \short Set vector of external halo node in this Mesh
2084  /// whose non-halo external counterpart is held on processor p.
2085  void set_external_halo_node_pt(const unsigned& p,
2087  {
2089  }
2090 
2091  /// Null out specified external halo node (used when deleting duplicates)
2092  void null_external_halo_node(const unsigned& p, Node* nod_pt);
2093 
2094  /// \short Consolidate external halo node storage by removing nulled out
2095  /// pointes in external halo and haloed schemes
2097 
2098  /// \short Total number of external haloed nodes in this Mesh
2100  {
2101  unsigned n=0;
2102  for (std::map<unsigned,Vector<Node*> >::iterator it=
2103  External_haloed_node_pt.begin();it!=External_haloed_node_pt.end();it++)
2104  {
2105  n+=it->second.size();
2106  }
2107  return n;
2108  }
2109 
2110  /// \short Number of external haloed nodes in this Mesh
2111  /// whose halo (external) counterpart is held on processor p.
2112  unsigned nexternal_haloed_node(const unsigned& p)
2113  {
2114  // Memory saving version of: return External_haloed_node_pt[p].size();
2115  std::map<unsigned, Vector<Node*> >::iterator it=
2116  External_haloed_node_pt.find(p);
2117  if (it==External_haloed_node_pt.end())
2118  {
2119  return 0;
2120  }
2121  return (*it).second.size();
2122  }
2123 
2124  /// \short Access fct to the j-th external haloed node in this Mesh
2125  /// whose halo external counterpart is held on processor p.
2126  Node* &external_haloed_node_pt(const unsigned &p, const unsigned &j)
2127  {
2128  return External_haloed_node_pt[p][j];
2129  }
2130 
2131  /// \short Add external haloed node whose halo (external) counterpart
2132  /// is held on processor p to the storage scheme for haloed nodes.
2133  unsigned add_external_haloed_node_pt(const unsigned& p, Node*& nod_pt);
2134 
2135  /// \short Access fct to vector of external haloed node in this Mesh
2136  /// whose halo external counterpart is held on processor p. (read only)
2138  {
2139  std::map<unsigned, Vector<Node*> >::iterator it=
2140  External_haloed_node_pt.find(p);
2141  if (it==External_haloed_node_pt.end())
2142  {
2143  Vector<Node*> tmp;
2144  return tmp;
2145  }
2146  return (*it).second;
2147  }
2148 
2149  /// \short Set vector of external haloed node in this Mesh
2150  /// whose halo external counterpart is held on processor p.
2151  void set_external_haloed_node_pt(const unsigned& p,
2153  {
2155  }
2156 
2157  /// \short Return the set of processors that hold external halo nodes. This is
2158  /// required to avoid having to pass a communicator into the node_update
2159  /// functions for Algebraic-based and MacroElement-based Meshes
2160  std::set<int> external_halo_proc()
2161  {
2162  std::set<int> procs;
2163  for (std::map<unsigned,Vector<Node*> >::iterator it=
2164  External_halo_node_pt.begin();it!=External_halo_node_pt.end();it++)
2165  {
2166  procs.insert((*it).first);
2167  }
2168  return procs;
2169  }
2170 
2171  /// \short Creates the shared boundaries, only used in unstructured meshes
2172  /// In this case with the "TriangleMesh" class
2174  const Vector<unsigned> &element_domain,
2176  &backed_up_el_pt,
2178  &backed_up_f_el_pt,
2179  std::map<Data*,std::set<unsigned> >
2180  &processors_associated_with_data,
2181  const bool&
2182  overrule_keep_as_halo_element_status)
2183  {
2184  std::ostringstream error_stream;
2185  error_stream << "Empty default create_shared_boundaries() method"
2186  << "called.\n";
2187  error_stream << "This should be overloaded in a specific "
2188  << "TriangleMeshBase\n";
2189  throw OomphLibError(error_stream.str(),
2190  "Mesh::create_shared_boundaries()",
2191  OOMPH_EXCEPTION_LOCATION);
2192  }
2193 
2194  // Check if necessary to add the element as haloed or if it has been
2195  // previously added to the haloed scheme
2196  virtual unsigned try_to_add_root_haloed_element_pt(const unsigned& p,
2197  GeneralisedElement*& el_pt)
2198  {
2199  std::ostringstream error_stream;
2200  error_stream << "Empty default try_to_add_root_haloed_element_pt() method"
2201  << "called.\n";
2202  error_stream << "This should be overloaded in a specific "
2203  << "TriangleMeshBase\n";
2204  throw OomphLibError(error_stream.str(),
2205  "Mesh::try_to_add_root_haloed_element_pt()",
2206  OOMPH_EXCEPTION_LOCATION);
2207  }
2208 
2209  // Check if necessary to add the node as haloed or if it has been
2210  // previously added to the haloed scheme
2211  virtual unsigned try_to_add_haloed_node_pt(const unsigned& p, Node*& nod_pt)
2212  {
2213  std::ostringstream error_stream;
2214  error_stream << "Empty default try_to_add_haloed_node_pt() method"
2215  << "called.\n";
2216  error_stream << "This should be overloaded in a specific "
2217  << "TriangleMeshBase\n";
2218  throw OomphLibError(error_stream.str(),
2219  "Mesh::try_to_add_haloed_node_pt()",
2220  OOMPH_EXCEPTION_LOCATION);
2221  }
2222 
2223 #endif
2224 
2225  /// \short Wipe the storage for all externally-based elements
2227 
2228 
2229 
2230 
2231 };
2232 
2233 
2234 
2235 //////////////////////////////////////////////////////////////////////////
2236 //////////////////////////////////////////////////////////////////////////
2237 //////////////////////////////////////////////////////////////////////////
2238 
2239 class SolidICProblem;
2240 
2241 //========================================================================
2242 /// \short General SolidMesh class.
2243 ///
2244 /// Solid meshes contain SolidFiniteElements which contain
2245 /// SolidNodes. This class simply overloads the appropriate access
2246 /// functions from the underlying Mesh class.
2247 //
2248 // Needs to be derived from Mesh with virtual so that
2249 // solid meshes can be derived from general meshes, without
2250 // multiple copies of Mesh objects
2251 //========================================================================
2252 class SolidMesh : public virtual Mesh
2253 {
2254  public:
2255 
2256 
2257  /// \short Default constructor
2259 
2260  /// Broken copy constructor
2261  SolidMesh(const SolidMesh& dummy)
2262  {
2263  BrokenCopy::broken_copy("SolidMesh");
2264  }
2265 
2266  /// Broken assignment operator
2267  void operator=(const SolidMesh&)
2268  {
2269  BrokenCopy::broken_assign("SolidMesh");
2270  }
2271 
2272 
2273  /// \short Constructor builds combined mesh from the meshes specified.
2274  /// Note: This simply merges the meshes' elements and nodes (ignoring
2275  /// duplicates; no boundary information etc. is created).
2276  SolidMesh(const Vector<SolidMesh*>& sub_mesh_pt)
2277  {
2278 #ifdef OOMPH_HAS_MPI
2279  // Mesh hasn't been distributed: Null out pointer to communicator
2280  Comm_pt=0;
2281 #endif
2282 
2283  unsigned n=sub_mesh_pt.size();
2284  Vector<Mesh*> sub_mesh_mesh_pt(n);
2285  for (unsigned i=0;i<n;i++)
2286  {
2287  sub_mesh_mesh_pt[i]=static_cast<Mesh*>(sub_mesh_pt[i]);
2288  }
2289  merge_meshes(sub_mesh_mesh_pt);
2290  }
2291 
2292  /// Return a pointer to the n-th global SolidNode
2293  //Can safely cast the nodes to SolidNodes
2294  SolidNode* node_pt(const unsigned long &n)
2295  {
2296 #ifdef PARANOID
2297  if(!dynamic_cast<SolidNode*>(Node_pt[n]))
2298  {
2299  std::ostringstream error_stream;
2300  error_stream << "Error: Node " << n << "is a "
2301  << typeid(Node_pt[n]).name()
2302  << ", not an SolidNode" << std::endl;
2303  throw OomphLibError(error_stream.str(),
2304  OOMPH_CURRENT_FUNCTION,
2305  OOMPH_EXCEPTION_LOCATION);
2306  }
2307 #endif
2308  //Return a static cast to the Node_pt
2309  return (static_cast<SolidNode*>(Node_pt[n]));
2310  }
2311 
2312  /// Return n-th SolidNodes on b-th boundary
2313  SolidNode* boundary_node_pt(const unsigned &b,
2314  const unsigned &n)
2315  {
2316 #ifdef PARANOID
2317  if(!dynamic_cast<SolidNode*>(Mesh::boundary_node_pt(b,n)))
2318  {
2319  std::ostringstream error_stream;
2320  error_stream
2321  << "Error: Node " << n << "is a "
2322  << typeid(Mesh::boundary_node_pt(b,n)).name()
2323  << ", not an SolidNode" << std::endl;
2324  throw OomphLibError(error_stream.str(),
2325  OOMPH_CURRENT_FUNCTION,
2326  OOMPH_EXCEPTION_LOCATION);
2327  }
2328 #endif
2329  return static_cast<SolidNode*>(Mesh::boundary_node_pt(b,n));
2330  }
2331 
2332  /// \short Return the n-th local SolidNode in elemnet e.
2333  ///This is required to cast the nodes in a solid mesh to be
2334  ///SolidNodes and therefore allow access to the extra SolidNode data
2335  SolidNode* element_node_pt(const unsigned long &e,
2336  const unsigned &n)
2337  {
2338 #ifdef PARANOID
2339  // Try to cast to FiniteElement
2340  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Element_pt[e]);
2341  if (el_pt==0)
2342  {
2343  //Error
2344  throw OomphLibError("Failed cast to FiniteElement* ",
2345  OOMPH_CURRENT_FUNCTION,
2346  OOMPH_EXCEPTION_LOCATION);
2347  }
2348  if(!dynamic_cast<SolidNode*>(el_pt->node_pt(n)))
2349  {
2350  std::ostringstream error_message;
2351  Node *np = el_pt->node_pt(n);
2352  error_message << "Error: Node " << n << " of element " << e
2353  << "is a " << typeid(*np).name()
2354  << ", not an SolidNode" << std::endl;
2355 
2356  throw OomphLibError(error_message.str(),
2357  OOMPH_CURRENT_FUNCTION,
2358  OOMPH_EXCEPTION_LOCATION);
2359  }
2360 #endif
2361  //Return a cast to an SolidNode
2362  return(static_cast<SolidNode*>(
2363  dynamic_cast<FiniteElement*>(Element_pt[e])->node_pt(n)));
2364  }
2365 
2366  /// \short Make the current configuration the undeformed one by
2367  /// setting the nodal Lagrangian coordinates to their current
2368  /// Eulerian ones
2370 
2371 
2372  /// \short Scale all nodal coordinates by given factor and re-assign the
2373  /// Lagrangian coordinates
2374  void scale_mesh(const double& factor)
2375  {
2376  Mesh::scale_mesh(factor);
2377 
2378  //Re-assign the Lagrangian coordinates
2380  }
2381  /// \short Static problem that can be used to assign initial conditions
2382  /// on a given solid mesh (need to define this as a static problem
2383  /// somewhere because deleting the problem would wipe out the mesh too!)
2385 
2386 
2387 };
2388 
2389 
2390 
2391 
2392 ///////////////////////////////////////////////////////////////////////////
2393 ///////////////////////////////////////////////////////////////////////////
2394 ///////////////////////////////////////////////////////////////////////////
2395 
2396 
2397 
2398 //===================================================
2399 /// Edge class
2400 //===================================================
2401  class Edge
2402  {
2403 
2404  public:
2405 
2406  /// Constructor: Pass in the two vertex nodes
2408  {
2409  if (node1_pt==node2_pt)
2410  {
2411 #ifdef PARANOID
2412  std::ostringstream error_stream;
2413  error_stream << "Edge cannot have two identical vertex nodes\n";
2414  throw OomphLibError(
2415  error_stream.str(),
2416  OOMPH_CURRENT_FUNCTION,
2417  OOMPH_EXCEPTION_LOCATION);
2418 #endif
2419  }
2420 
2421 
2422  // Sort lexicographically based on pointer address of nodes
2423  if (node1_pt>node2_pt)
2424  {
2427  }
2428  else
2429  {
2432  }
2433  }
2434 
2435 
2436  /// Access to the first vertex node
2437  Node* node1_pt() const {return Node1_pt;}
2438 
2439  /// Access to the second vertex node
2440  Node* node2_pt() const {return Node2_pt;}
2441 
2442  /// Comparison operator
2443  bool operator==(const Edge& other) const
2444  {
2445  if ((Node1_pt==other.node1_pt())&&
2446  (Node2_pt==other.node2_pt()))
2447 
2448  {
2449  return true;
2450  }
2451  else
2452  {
2453  return false;
2454  }
2455  }
2456 
2457 
2458 
2459  /// Less-than operator
2460  bool operator<(const Edge& other) const
2461  {
2462  if (Node1_pt<other.node1_pt())
2463  {
2464  return true;
2465  }
2466  else if (Node1_pt==other.node1_pt())
2467  {
2468  if (Node2_pt<other.node2_pt())
2469  {
2470  return true;
2471  }
2472  else
2473  {
2474  return false;
2475  }
2476  }
2477  else
2478  {
2479  return false;
2480  }
2481  }
2482 
2483  /// \short Test whether the Edge lies on a boundary. Relatively simple
2484  /// test, based on both vertices lying on (some) boundary.
2485  bool is_on_boundary() const
2486  {
2487  return (Node1_pt->is_on_boundary() && Node2_pt->is_on_boundary() );
2488  }
2489 
2490 
2491  /// \short Test whether the Edge is a boundary edge, i.e. does it
2492  /// connnect two boundary nodes?
2493  bool is_boundary_edge() const
2494  {
2495  return ((dynamic_cast<BoundaryNodeBase*>(Node1_pt)!=0)&&
2496  (dynamic_cast<BoundaryNodeBase*>(Node2_pt)!=0));
2497  }
2498 
2499  private:
2500 
2501  /// First vertex node
2503 
2504  /// Second vertex node
2506  };
2507 
2508 
2509 
2510 ///////////////////////////////////////////////////////////////////////
2511 ///////////////////////////////////////////////////////////////////////
2512 ///////////////////////////////////////////////////////////////////////
2513 
2514 
2515 
2516 //=================================================================
2517 /// Namespace with helper function to check element type in mesh
2518 /// constructors (say).
2519 //=================================================================
2520 namespace MeshChecker
2521 {
2522 
2523 
2524 //=================================================================
2525 /// \short Helper function to assert that finite element of type ELEMENT
2526 /// can be cast to base class of type GEOM_ELEMENT_BASE -- useful
2527 /// to avoid confusion if a mesh that was written for a specific
2528 /// element type (e.g. a QElement) is used with another one (e.g.
2529 /// a TElement. First argument specifies the required spatial dimension
2530 /// of the element (i.e. the number of local coordinates). The optional
2531 /// second argument specifies the required nnode_1d (i.e. the number
2532 /// of nodes along a 1D element edge). Can be omitted if the mesh
2533 /// can handle any number in which case this test is skipped.
2534 //=================================================================
2535  template<class GEOM_ELEMENT_BASE, class ELEMENT>
2536  void assert_geometric_element(const unsigned& dim,
2537  const unsigned& nnode_1d=0)
2538  {
2539  // Only do tests in paranoia mode
2540 #ifndef PARANOID
2541  return;
2542 #endif
2543 
2544  // Instantiate element
2545  ELEMENT* el_pt=new ELEMENT;
2546 
2547  // Can we cast to required geometric element base type
2548  if (dynamic_cast<GEOM_ELEMENT_BASE*>(el_pt)==0)
2549  {
2550  std::stringstream error_message;
2551  error_message
2552  << "You have specified an illegal element type! Element is of type \n\n"
2553  << typeid(el_pt).name()
2554  << "\n\nand cannot be cast to type \n\n "
2555  << typeid(GEOM_ELEMENT_BASE).name()
2556  << "\n\n";
2557  throw OomphLibError(error_message.str(),
2558  OOMPH_CURRENT_FUNCTION,
2559  OOMPH_EXCEPTION_LOCATION);
2560  }
2561 
2562  // Does the dimension match?
2563  if (dim!=el_pt->dim())
2564  {
2565  std::stringstream error_message;
2566  error_message
2567  << "You have specified an illegal element type! Element is of type \n\n"
2568  << typeid(el_pt).name()
2569  << "\n\nand has dimension = " << el_pt->dim()
2570  << " but we need dim = " << dim << std::endl;
2571  throw OomphLibError(error_message.str(),
2572  OOMPH_CURRENT_FUNCTION,
2573  OOMPH_EXCEPTION_LOCATION);
2574 
2575  }
2576 
2577  // Does nnode_1d match?
2578  if (nnode_1d!=0)
2579  {
2580  if (nnode_1d!=el_pt->nnode_1d())
2581  {
2582  std::stringstream error_message;
2583  error_message
2584  << "You have specified an illegal element type! Element is of type \n\n"
2585  << typeid(el_pt).name()
2586  << "\n\nand has nnode_1d = " << el_pt->nnode_1d()
2587  << " but we need nnode_1d = " << nnode_1d << std::endl;
2588  throw OomphLibError(error_message.str(),
2589  OOMPH_CURRENT_FUNCTION,
2590  OOMPH_EXCEPTION_LOCATION);
2591  }
2592  }
2593 
2594  // Clean up
2595  delete el_pt;
2596 
2597  }
2598 
2599 }
2600 
2601 
2602 
2603 ///////////////////////////////////////////////////////////////////////
2604 ///////////////////////////////////////////////////////////////////////
2605 ///////////////////////////////////////////////////////////////////////
2606 
2607 
2608 //=================paraview_helper=================================
2609 /// Namespace for paraview-style output helper functions
2610 //=================================================================
2611 namespace ParaviewHelper
2612 {
2613 
2614  /// Write the pvd file header
2615  extern void write_pvd_header(std::ofstream &pvd_file);
2616 
2617  /// \short Add name of output file and associated continuous time
2618  /// to pvd file.
2619  extern void write_pvd_information(std::ofstream &pvd_file,
2620  const std::string& output_filename,
2621  const double& time);
2622 
2623  /// Write the pvd file footer
2624  extern void write_pvd_footer(std::ofstream &pvd_file);
2625 
2626 }
2627 
2628 ////////////////////////////////////////////////////////////////
2629 ////////////////////////////////////////////////////////////////
2630 ////////////////////////////////////////////////////////////////
2631 
2632 namespace NodeOrdering
2633 {
2634  /// Function for ordering nodes. Return true if first node's position is
2635  /// "before" second nodes. Dimension 0 checked first, then... until they
2636  /// are different (by more than tol=1e-10). If they are both in exactly
2637  /// the same place an error is thrown.
2638  inline bool node_global_position_comparison(Node* nd1_pt, Node* nd2_pt)
2639  {
2640  unsigned ndim = nd1_pt->ndim();
2641 
2642  unsigned j;
2643  for(j=0; j<ndim; j++)
2644  {
2645  if(std::abs(nd1_pt->x(j) - nd2_pt->x(j)) > 1e-10)
2646  {
2647  if(nd1_pt->x(j) < nd2_pt->x(j))
2648  {
2649  return true;
2650  }
2651  else
2652  {
2653  return false;
2654  }
2655  }
2656  // otherwise they are the same in this dimension, check next one.
2657  }
2658 
2659  std::string err = "Nodes are at the same point to ~ 1e-10!";
2660  err += " difference is " +
2661  StringConversion::to_string(std::abs(nd1_pt->x(j) - nd2_pt->x(j)));
2662  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
2663  OOMPH_CURRENT_FUNCTION);
2664  }
2665 }
2666 
2667 
2668 }
2669 
2670 #endif
unsigned nshared_node()
Total number of shared nodes in this Mesh.
Definition: mesh.h:1756
static bool Suppress_warning_about_empty_mesh_level_time_stepper_function
Boolean used to control warning about empty mesh level timestepper function.
Definition: mesh.h:246
void output_boundaries(const std::string &output_filename)
Definition: mesh.h:970
void assign_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers in all elements If the boolean argument is true then also store poi...
Definition: mesh.cc:702
A Generalised Element class.
Definition: elements.h:76
void set_communicator_pt(OomphCommunicator *comm_pt)
Definition: mesh.h:1332
Node * haloed_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th haloed node in this Mesh whose halo counterpart is held on processor p...
Definition: mesh.h:1712
void doc_boundary_coordinates(const unsigned &b, std::ofstream &the_file)
Output boundary coordinates on boundary b – template argument specifies the bulk element type (needed...
Definition: mesh.h:316
void add_shared_node_pt(const unsigned &p, Node *&nod_pt)
Add shared node whose counterpart is held on processor p to the storage scheme for shared nodes...
Definition: mesh.h:1828
GeneralisedElement *& root_haloed_element_pt(const unsigned &p, const unsigned &e)
Access fct to the e-th root haloed element in this Mesh whose non-halo counterpart is held on process...
Definition: mesh.h:1666
unsigned nexternal_halo_node(const unsigned &p)
Number of external halo nodes in this Mesh whose non-halo (external) counterpart is held on processor...
Definition: mesh.h:2041
bool operator==(const Edge &other) const
Comparison operator.
Definition: mesh.h:2443
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
unsigned nexternal_haloed_node(const unsigned &p)
Number of external haloed nodes in this Mesh whose halo (external) counterpart is held on processor p...
Definition: mesh.h:2112
void disable_resizing_of_halo_nodes()
Function to suppress resizing of halo nodes – optmisation but call it at your own risk! ...
Definition: mesh.h:1730
std::map< unsigned, Vector< Node * > > Shared_node_pt
Definition: mesh.h:127
void write_pvd_footer(std::ofstream &pvd_file)
Write the pvd file footer.
Definition: mesh.cc:9281
Vector< GeneralisedElement * > halo_element_pt(const unsigned &p)
Return vector of halo elements in this Mesh whose non-halo counterpart is held on processor p...
Definition: mesh.h:1445
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:194
std::map< unsigned, Vector< GeneralisedElement * > > Root_halo_element_pt
Map of vectors holding the pointers to the root halo elements.
Definition: mesh.h:111
void output_boundaries_coordinates(std::ostream &outfile)
Definition: mesh.cc:1025
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:497
unsigned nexternal_halo_element()
Total number of external halo elements in this Mesh.
Definition: mesh.h:1913
std::map< unsigned, Vector< Node * > > Haloed_node_pt
Map of vectors holding the pointers to the haloed nodes.
Definition: mesh.h:120
SolidMesh(const SolidMesh &dummy)
Broken copy constructor.
Definition: mesh.h:2261
SolidNode * element_node_pt(const unsigned long &e, const unsigned &n)
Return the n-th local SolidNode in elemnet e. This is required to cast the nodes in a solid mesh to b...
Definition: mesh.h:2335
void doc_shared_nodes()
Doc shared nodes.
Definition: mesh.h:1768
std::map< unsigned, Vector< GeneralisedElement * > > External_halo_element_pt
Map of vectors holding the pointers to the external halo elements.
Definition: mesh.h:137
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:246
void add_external_halo_node_pt(const unsigned &p, Node *&nod_pt)
Add external halo node whose non-halo (external) counterpart is held on processor p to the storage sc...
Definition: mesh.h:2055
void resize_halo_nodes()
Helper function that resizes halo nodes to the same size as their non-halo counterparts if required...
Definition: mesh.cc:4126
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors. ...
Definition: mesh.h:85
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
SolidMesh()
Default constructor.
Definition: mesh.h:2258
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt)
Output a given Vector function at f(n_plot) points in each element.
Definition: mesh.cc:1953
std::vector< bool > Boundary_coordinate_exists
Vector of boolean data that indicates whether the boundary coordinates have been set for the boundary...
Definition: mesh.h:201
unsigned add_external_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add external haloed element whose non-halo counterpart is held on processor p to the storage scheme f...
Definition: mesh.cc:9107
Vector< GeneralisedElement * > & element_pt()
Return reference to the Vector of elements.
Definition: mesh.h:473
virtual ~Mesh()
Virtual Destructor to clean up all memory.
Definition: mesh.cc:588
void prune_halo_elements_and_nodes(Vector< GeneralisedElement * > &deleted_element_pt, const bool &report_stats=false)
(Irreversibly) prune halo(ed) elements and nodes, usually after another round of refinement, to get rid of excessively wide halo layers. Note that the current mesh will be now regarded as the base mesh and no unrefinement relative to it will be possible once this function has been called.
Definition: mesh.h:1372
Information for documentation of results: Directory and file number to enable output in the form RESL...
Vector< Node * > external_halo_node_pt(const unsigned &p)
Access fct to vector of external halo node in this Mesh whose non-halo external counterpart is held o...
Definition: mesh.h:2071
Mesh()
Default constructor.
Definition: mesh.h:249
cstr elem_len * i
Definition: cfortran.h:607
Mesh(const Vector< Mesh * > &sub_mesh_pt)
Constructor builds combined mesh from the meshes specified. Note: This simply merges the meshes' elem...
Definition: mesh.h:270
void add_root_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add root haloed element whose non-halo counterpart is held on processor p to the storage scheme for h...
Definition: mesh.h:1677
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
void add_halo_node_pt(const unsigned &p, Node *&nod_pt)
Add halo node whose non-halo counterpart is held on processor p to the storage scheme for halo nodes...
Definition: mesh.h:1608
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
void shift_time_values()
Shift time-dependent data along for next timestep: Deal with nodal Data/positions and the element's i...
Definition: mesh.cc:2076
std::map< unsigned, Vector< Node * > > External_haloed_node_pt
Map of vectors holding the pointers to the external haloed nodes.
Definition: mesh.h:149
void output_paraview(std::ofstream &file_out, const unsigned &nplot) const
Output in paraview format into specified file. Breaks up each element into sub-elements for plotting ...
Definition: mesh.cc:1184
GeneralisedElement *& external_halo_element_pt(const unsigned &p, const unsigned &e)
Access fct to the e-th external halo element in this Mesh whose non-halo counterpart is held on proce...
Definition: mesh.h:1941
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1305
void get_shared_node_pt(const unsigned &p, Vector< Node * > &shared_node_pt)
Get vector of pointers to shared nodes with processor p. Required for faster search in Missing_master...
Definition: mesh.h:1814
void add_root_halo_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add root halo element whose non-halo counterpart is held on processor p to this Mesh.
Definition: mesh.h:1574
Node * get_some_non_boundary_node() const
Find a node not on any boundary in mesh_pt (useful for pinning a single node in a purely Neumann prob...
Definition: mesh.h:832
The Problem class.
Definition: problem.h:152
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
unsigned nexternal_halo_node()
Total number of external halo nodes in this Mesh.
Definition: mesh.h:1997
virtual void distribute(OomphCommunicator *comm_pt, const Vector< unsigned > &element_domain, Vector< GeneralisedElement * > &deleted_element_pt, DocInfo &doc_info, const bool &report_stats, const bool &overrule_keep_as_halo_element_status)
Distribute the problem and doc; make this virtual to allow overloading for particular meshes where fu...
Definition: mesh.cc:4640
bool Output_halo_elements
Bool for output of halo elements.
Definition: mesh.h:1725
A general Finite Element class.
Definition: elements.h:1271
void enable_output_of_halo_elements()
Function to enable halo element output.
Definition: mesh.h:1750
bool is_boundary_edge() const
Test whether the Edge is a boundary edge, i.e. does it connnect two boundary nodes?
Definition: mesh.h:2493
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary, so return false. This will be overloaded by BoundaryNodes.
Definition: nodes.h:1293
void get_halo_node_stats(double &av_number, unsigned &max_number, unsigned &min_number)
Get halo node stats for this distributed mesh: Average/max/min number of halo nodes over all processo...
Definition: mesh.cc:4534
void operator=(const Mesh &)
Broken assignment operator.
Definition: mesh.h:408
virtual void set_mesh_level_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Function that can be used to set any additional timestepper data stored at the Mesh (as opposed to no...
Definition: mesh.cc:2149
Node * node1_pt() const
Access to the first vertex node.
Definition: mesh.h:2437
void unset_keep_all_elements_as_halos()
Calll this function to unset the flag that keeps all elements in the mesh as halo elements...
Definition: mesh.h:1339
virtual void scale_mesh(const double &factor)
Scale all nodal coordinates by given factor. Virtual so it can be overloaded in SolidMesh class where...
Definition: mesh.h:386
void output_fct_paraview(std::ofstream &file_out, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
Output in paraview format into specified file. Breaks up each element into sub-elements for plotting ...
Definition: mesh.cc:1482
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
unsigned self_test()
Self-test: Check elements and nodes. Return 0 for OK.
Definition: mesh.cc:715
Vector< GeneralisedElement * > haloed_element_pt(const unsigned &p)
Return vector of haloed elements in this Mesh whose haloing counterpart is held on processor p...
Definition: mesh.h:1484
void operator=(const SolidMesh &)
Broken assignment operator.
Definition: mesh.h:2267
GeneralisedElement *& root_halo_element_pt(const unsigned &p, const unsigned &e)
Access fct to the e-th root halo element in this Mesh whose non-halo counterpart is held on processor...
Definition: mesh.h:1565
OomphInfo oomph_info
void set_elemental_internal_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set the timestepper associated with the internal data stored within elements in the meah...
Definition: mesh.cc:2279
virtual void dump(std::ofstream &dump_file, const bool &use_old_ordering=true) const
Dump the data in the mesh into a file for restart.
Definition: mesh.cc:1044
unsigned nroot_haloed_element()
Total number of root haloed elements in this Mesh.
Definition: mesh.h:1623
void output_external_haloed_elements(std::ostream &outfile, const unsigned &n_plot=5)
Output all external haloed elements.
Definition: mesh.h:1883
unsigned nhalo_node(const unsigned &p)
Number of halo nodes in this Mesh whose non-halo counterpart is held on processor p...
Definition: mesh.h:1594
void build_face_mesh(const unsigned &b, Mesh *const &face_mesh_pt)
Constuct a Mesh of FACE_ELEMENTs along the b-th boundary of the mesh (which contains elements of type...
Definition: mesh.h:639
virtual void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
Plot error when compared against a given time-depdendent exact solution. Also returns the norm of the...
Definition: mesh.h:1244
e
Definition: cfortran.h:575
void set_boundary_number_in_bulk_mesh(const unsigned &b)
Set function for the boundary number in bulk mesh.
Definition: elements.h:4236
unsigned nnon_halo_element()
Total number of non-halo elements in this mesh (Costly call computes result on the fly) ...
Definition: mesh.h:1522
bool Lookup_for_elements_next_boundary_is_setup
Definition: mesh.h:98
Node * node_pt(const unsigned long &n) const
Return pointer to global node n (const version)
Definition: mesh.h:459
void get_efficiency_of_mesh_distribution(double &av_efficiency, double &max_efficiency, double &min_efficiency)
Get efficiency of mesh distribution: In an ideal distribution without halo overhead, each processor would only hold its own elements. Efficieny per processor = (number of non-halo elements)/ (total number of elements).
Definition: mesh.cc:6070
void check_inverted_elements(bool &mesh_has_inverted_elements)
Check for inverted elements and report outcome.
Definition: mesh.h:727
void convert_to_boundary_node(Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
A function that upgrades an ordinary node to a boundary node We shouldn't ever really use this...
Definition: mesh.cc:2334
unsigned nexternal_haloed_element(const unsigned &p)
Number of external haloed elements in this Mesh whose non-halo counterpart is held on processor p...
Definition: mesh.h:1971
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:809
void copy_boundary_node_data_from_nodes()
Definition: mesh.h:526
void get_haloed_node_stats(double &av_number, unsigned &max_number, unsigned &min_number)
Get haloed node stats for this distributed mesh: Average/max/min number of haloed nodes over all proc...
Definition: mesh.cc:4589
void write_pvd_information(std::ofstream &pvd_file, const std::string &output_filename, const double &time)
Add name of output file and associated continuous time to pvd file.
Definition: mesh.cc:9260
void remove_boundary_node(const unsigned &b, Node *const &node_pt)
Definition: mesh.cc:222
virtual void compute_norm(double &norm)
Compute norm of solution – broken virtual can be overloaded by element writer to implement whatever n...
Definition: elements.h:1119
void output(const std::string &output_filename, const unsigned &n_plot)
Output at f(n_plot) points in each element.
Definition: mesh.h:947
virtual void classify_halo_and_haloed_nodes(const bool &report_stats=false)
Definition: mesh.h:1418
void get_external_halo_node_pt(Vector< Node * > &external_halo_node_pt)
Get vector of pointers to all external halo nodes.
Definition: mesh.h:2010
void output_boundary_coordinates(const unsigned &b, const std::string &output_filename)
Definition: mesh.h:994
void output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:1778
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
void assign_initial_values_impulsive()
Assign initial values for an impulsive start.
Definition: mesh.cc:2040
Vector< Node * > external_haloed_node_pt(const unsigned &p)
Access fct to vector of external haloed node in this Mesh whose halo external counterpart is held on ...
Definition: mesh.h:2137
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:602
void get_all_halo_data(std::map< unsigned, double * > &map_of_halo_data)
Get all the halo data stored in the mesh and add pointers to the data to the map, indexed by global e...
Definition: mesh.cc:4438
unsigned nodal_dimension() const
Return number of nodal dimension in mesh.
Definition: mesh.cc:8699
void flush_element_storage()
Flush storage for elements (only) by emptying the vectors that store the pointers to them...
Definition: mesh.h:443
SolidNode * node_pt(const unsigned long &n)
Return a pointer to the n-th global SolidNode.
Definition: mesh.h:2294
void scale_mesh(const double &factor)
Scale all nodal coordinates by given factor and re-assign the Lagrangian coordinates.
Definition: mesh.h:2374
void dump(const std::string &dump_file_name, const bool &use_old_ordering=true) const
Dump the data in the mesh into a file for restart.
Definition: mesh.h:889
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Return a pointer to set of mesh boundaries that this node occupies; this will be overloaded by Bounda...
Definition: nodes.h:1287
unsigned nexternal_halo_element(const unsigned &p)
Number of external halo elements in this Mesh whose non-halo counterpart is held on processor p...
Definition: mesh.h:1927
Node * node2_pt() const
Access to the second vertex node.
Definition: mesh.h:2440
GeneralisedElement * element_pt(const unsigned long &e) const
Return pointer to element e (const version)
Definition: mesh.h:466
void calculate_predictions()
Calculate predictions for all Data and positions associated with the mesh, usually used in adaptive t...
Definition: mesh.cc:2115
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
unsigned add_external_haloed_node_pt(const unsigned &p, Node *&nod_pt)
Add external haloed node whose halo (external) counterpart is held on processor p to the storage sche...
Definition: mesh.cc:9148
OomphCommunicator * communicator_pt() const
Definition: mesh.h:1317
double size() const
Definition: elements.cc:4121
void merge_meshes(const Vector< Mesh * > &sub_mesh_pt)
Merge meshes. Note: This simply merges the meshes' elements and nodes (ignoring duplicates; no bounda...
Definition: mesh.cc:69
bool Keep_all_elements_as_halos
bool to indicate whether to keep all elements in a mesh as halos or not
Definition: mesh.h:152
void check_inverted_elements(bool &mesh_has_inverted_elements, std::ofstream &inverted_element_file)
Check for inverted elements and report outcome in boolean variable. This visits all elements at their...
Definition: mesh.cc:802
void null_external_halo_node(const unsigned &p, Node *nod_pt)
Null out specified external halo node (used when deleting duplicates)
Definition: mesh.cc:8226
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
void set_nodal_and_elemental_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set the timestepper associated with all nodal and elemental data stored in the mesh.
Definition: mesh.h:1032
bool operator<(const Edge &other) const
Less-than operator.
Definition: mesh.h:2460
unsigned nhalo_node()
Total number of halo nodes in this Mesh.
Definition: mesh.h:1581
void set_keep_all_elements_as_halos()
Call this function to keep all the elements as halo elements.
Definition: mesh.h:1335
virtual bool is_a_copy() const
Return a boolean to indicate whether the Data objact contains any copied values. A base Data object c...
Definition: nodes.h:255
std::map< unsigned, Vector< GeneralisedElement * > > Root_haloed_element_pt
Map of vectors holding the pointers to the root haloed elements.
Definition: mesh.h:114
virtual unsigned try_to_add_root_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Definition: mesh.h:2196
std::map< unsigned, Vector< GeneralisedElement * > > External_haloed_element_pt
Map of vectors holding the pointers to the external haloed elements.
Definition: mesh.h:140
void flush_node_storage()
Flush storage for nodes (only) by emptying the vectors that store the pointers to them...
Definition: mesh.h:450
Edge(Node *node1_pt, Node *node2_pt)
Constructor: Pass in the two vertex nodes.
Definition: mesh.h:2407
Mesh(const Mesh &dummy)
Broken copy constructor.
Definition: mesh.h:402
Vector< Vector< int > > Face_index_at_boundary
For the e-th finite element on boundary b, this is the index of the face that lies along that boundar...
Definition: mesh.h:106
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
unsigned nhaloed_node(const unsigned &p)
Number of haloed nodes in this Mesh whose haloed counterpart is held on processor p...
Definition: mesh.h:1698
Node *& external_halo_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th external halo node in this Mesh whose non-halo external counterpart is held on...
Definition: mesh.h:2064
void set_external_halo_node_pt(const unsigned &p, const Vector< Node * > &external_halo_node_pt)
Set vector of external halo node in this Mesh whose non-halo external counterpart is held on processo...
Definition: mesh.h:2085
unsigned nexternal_haloed_element()
Total number of external haloed elements in this Mesh.
Definition: mesh.h:1957
void set_halo(const unsigned &non_halo_proc_ID)
Label the element as halo and specify processor that holds non-halo counterpart.
Definition: elements.h:1132
static char t char * s
Definition: cfortran.h:572
double total_size()
Determine the sum of all "sizes" of the FiniteElements in the mesh (non-FiniteElements are ignored)...
Definition: mesh.h:693
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
virtual void classify_halo_and_haloed_nodes(DocInfo &doc_info, const bool &report_stats)
Classify the halo and haloed nodes in the mesh. Virtual so it can be overloaded to perform additional...
Definition: mesh.cc:3004
unsigned nexternal_haloed_node()
Total number of external haloed nodes in this Mesh.
Definition: mesh.h:2099
Tree * tree_pt()
Access function: Pointer to quadtree representation of this element.
Vector< Node * > prune_dead_nodes()
Prune nodes. Nodes that have been marked as obsolete are removed from the mesh (and its boundary-node...
Definition: mesh.cc:900
Node * Node1_pt
First vertex node.
Definition: mesh.h:2502
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
Definition: mesh.h:470
Node * shared_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th shared node in this Mesh who has a counterpart on processor p...
Definition: mesh.h:1805
void output_boundaries(std::ostream &outfile)
Output the nodes on the boundaries (into separate tecplot zones)
Definition: mesh.cc:1001
void describe_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the dofs of the Mesh. The ostream specifies the output stream to which the descr...
Definition: mesh.cc:648
void distribute(OomphCommunicator *comm_pt, const Vector< unsigned > &element_domain, Vector< GeneralisedElement * > &deleted_element_pt, const bool &report_stats=false)
Distribute the problem Add to vector of pointers to deleted elements.
Definition: mesh.h:1353
virtual void compute_norm(double &norm)
Compute norm of solution by summing contributions of compute_norm(...) for all constituent elements i...
Definition: mesh.h:1068
void set_nodal_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set the timestepper associated with the nodal data in the mesh.
Definition: mesh.cc:2261
bool boundary_coordinate_exists(const unsigned &i) const
Indicate whether the i-th boundary has an intrinsic coordinate.
Definition: mesh.h:565
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:814
virtual void setup_boundary_element_info(std::ostream &outfile)
Setup lookup schemes which establish whic elements are located next to mesh's boundaries. Doc in outfile (if it's open). (Empty virtual function – implement this for specific Mesh classes)
Definition: mesh.h:294
Vector< GeneralisedElement * > root_halo_element_pt(const unsigned &p)
Vector of pointers to root halo elements in this Mesh whose non-halo counterpart is held on processor...
Definition: mesh.h:1557
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void add_haloed_node_pt(const unsigned &p, Node *&nod_pt)
Add haloed node whose halo counterpart is held on processor p to the storage scheme for haloed nodes...
Definition: mesh.h:1719
void remove_null_pointers_from_external_halo_node_storage()
Consolidate external halo node storage by removing nulled out pointes in external halo and haloed sch...
Definition: mesh.cc:8246
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:197
Node *& external_haloed_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th external haloed node in this Mesh whose halo external counterpart is held on p...
Definition: mesh.h:2126
void check_halo_schemes(DocInfo &doc_info, double &max_permitted_error_for_halo_check)
Check halo and shared schemes on the mesh.
Definition: mesh.cc:6571
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the elements. The ostream specifies the output stream to which...
Definition: mesh.cc:683
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:852
virtual void reset_boundary_element_info(Vector< unsigned > &ntmp_boundary_elements, Vector< Vector< unsigned > > &ntmp_boundary_elements_in_region, Vector< FiniteElement * > &deleted_elements)
Virtual function to perform the reset boundary elements info rutines.
Definition: mesh.h:297
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
For the e-th finite element on boundary b, return int to indicate the face_index of the face adjacent...
Definition: mesh.h:870
unsigned nhaloed_node()
Total number of haloed nodes in this Mesh.
Definition: mesh.h:1684
std::map< unsigned, Vector< Node * > > Halo_node_pt
Map of vectors holding the pointers to the halo nodes.
Definition: mesh.h:117
Edge class.
Definition: mesh.h:2401
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:456
GeneralisedTimestepper used to store the arclength derivatives and pervious solutions required in con...
Vector< Vector< FiniteElement * > > Boundary_element_pt
Vector of Vector of pointers to elements on the boundaries: Boundary_element_pt(b,e)
Definition: mesh.h:102
Vector< GeneralisedElement * > root_haloed_element_pt(const unsigned &p)
Vector of pointers to root haloed elements in this Mesh whose non-halo counterpart is held on process...
Definition: mesh.h:1651
Node * boundary_node_pt(const unsigned &b, const unsigned &n) const
Return pointer to node n on boundary b.
Definition: mesh.h:501
void doc_mesh_distribution(DocInfo &doc_info)
Doc the mesh distribution, to be processed with tecplot macros.
Definition: mesh.cc:6143
bool Resize_halo_nodes_not_required
Set this to true to suppress resizing of halo nodes (at your own risk!)
Definition: mesh.h:155
unsigned nroot_haloed_element(const unsigned &p)
Number of root haloed elements in this Mesh whose non-halo counterpart is held on processor p...
Definition: mesh.h:1636
void output_external_haloed_elements(const unsigned &p, std::ostream &outfile, const unsigned &n_plot=5)
Output all external haloed elements with processor p.
Definition: mesh.h:1896
std::map< unsigned, Vector< Node * > > External_halo_node_pt
Map of vectors holding the pointers to the external halo nodes.
Definition: mesh.h:146
General SolidMesh class.
Definition: mesh.h:2252
void output_external_halo_elements(const unsigned &p, std::ostream &outfile, const unsigned &n_plot=5)
Output all external halo elements with processor p.
Definition: mesh.h:1866
void output_boundaries_coordinates(const std::string &output_filename)
Definition: mesh.h:1009
void set_halo(const unsigned &non_halo_proc_ID)
Label the node as halo and specify processor that holds non-halo counterpart.
Definition: nodes.h:481
void disable_doc()
Disable documentation.
unsigned long assign_global_eqn_numbers(Vector< double * > &Dof_pt)
Assign the global equation numbers in the Data stored at the nodes and also internal element Data...
Definition: mesh.cc:614
void output(const std::string &output_filename)
Output for all elements.
Definition: mesh.h:938
virtual void node_update(const bool &update_all_solid_nodes=false)
Update nodal positions in response to changes in the domain shape. Uses the FiniteElement::get_x(...) function for FiniteElements and doesn't do anything for other element types. If a MacroElement pointer has been set for a FiniteElement, the MacroElement representation is used to update the nodal positions; if not get_x(...) uses the FE interpolation and thus leaves the nodal positions unchanged. Virtual, so it can be overloaded by specific meshes, such as AlgebraicMeshes or SpineMeshes. Generally, this function updates the position of all nodes in response to changes in the boundary position. However, we ignore all SolidNodes since their position is computed as part of the solution – unless the bool flag is set to true. Such calls are typically made when the initial mesh is created and/or after a mesh has been refined repeatedly before the start of the computation.
Definition: mesh.cc:290
virtual void output_boundary_coordinates(const unsigned &b, std::ostream &outfile)
Definition: mesh.h:980
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
Definition: nodes.h:1568
unsigned nroot_halo_element()
Total number of root halo elements in this Mesh.
Definition: mesh.h:1535
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4277
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
bool is_on_boundary() const
Test whether the Edge lies on a boundary. Relatively simple test, based on both vertices lying on (so...
Definition: mesh.h:2485
Vector< Vector< Node * > > Boundary_node_pt
Vector of Vector of pointers to nodes on the boundaries: Boundary_node_pt(b,n). Note that this is pri...
Definition: mesh.h:94
void synchronise_shared_nodes(const bool &report_stats)
Synchronise shared node lookup schemes to cater for the the case where: (1) a certain node on the cur...
Definition: mesh.cc:2635
void set_consistent_pinned_values_for_continuation(ContinuationStorageScheme *const &continuation_stepper_pt)
Set consistent values for pinned data in continuation.
Definition: mesh.cc:2184
void max_and_min_element_size(double &max_size, double &min_size)
Determine max and min area for all FiniteElements in the mesh (non-FiniteElements are ignored) ...
Definition: mesh.h:676
unsigned elemental_dimension() const
Return number of elemental dimension in mesh.
Definition: mesh.cc:8559
void output_external_halo_elements(std::ostream &outfile, const unsigned &n_plot=5)
Output all external halo elements.
Definition: mesh.h:1853
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned check_for_repeated_nodes(const double &epsilon=1.0e-12)
Check for repeated nodes within a given spatial tolerance. Return (0/1) for (pass/fail).
Definition: mesh.h:737
bool does_pointer_correspond_to_mesh_data(double *const &parameter_pt)
Does the double pointer correspond to any mesh data.
Definition: mesh.cc:2218
virtual void setup_boundary_element_info()
Interface for function that is used to setup the boundary information (Empty virtual function – imple...
Definition: mesh.h:288
static SolidICProblem Solid_IC_problem
Static problem that can be used to assign initial conditions on a given solid mesh (need to define th...
Definition: mesh.h:2384
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition: elements.cc:4480
void(FiniteElement::* SteadyExactSolutionFctPt)(const Vector< double > &x, Vector< double > &soln)
Typedef for function pointer to function that computes steady exact solution.
Definition: mesh.h:236
void add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
Definition: mesh.h:605
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
virtual void reorder_nodes(const bool &use_old_ordering=true)
Re-order nodes in the order in which they appear in elements – can be overloaded for more efficient r...
Definition: mesh.cc:483
bool node_global_position_comparison(Node *nd1_pt, Node *nd2_pt)
Definition: mesh.h:2638
unsigned nroot_halo_element(const unsigned &p)
Number of root halo elements in this Mesh whose non-halo counterpart is held on processor p...
Definition: mesh.h:1549
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
unsigned nshared_node(const unsigned &p)
Number of shared nodes in this Mesh who have a counterpart on processor p.
Definition: mesh.h:1791
virtual void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
Plot error when compared against a given time-depdendent exact solution. Also returns the norm of the...
Definition: mesh.h:1179
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:806
Node * Node2_pt
Second vertex node.
Definition: mesh.h:2505
virtual void get_node_reordering(Vector< Node * > &reordering, const bool &use_old_ordering=true) const
Get a reordering of the nodes in the order in which they appear in elements – can be overloaded for m...
Definition: mesh.cc:500
void disable_output_of_halo_elements()
Function to disable halo element output.
Definition: mesh.h:1744
void assert_geometric_element(const unsigned &dim, const unsigned &nnode_1d=0)
Helper function to assert that finite element of type ELEMENT can be cast to base class of type GEOM_...
Definition: mesh.h:2536
virtual void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Plot error when compared against a given time-depdendent exact solution. Also returns the norm of the...
Definition: mesh.h:1138
void flush_element_and_node_storage()
Flush storage for elements and nodes by emptying the vectors that store the pointers to them...
Definition: mesh.h:427
void enable_resizing_of_halo_nodes()
Function to (re-)enable resizing of halo nodes – this returns things to the default behaviour...
Definition: mesh.h:1738
virtual unsigned try_to_add_haloed_node_pt(const unsigned &p, Node *&nod_pt)
Definition: mesh.h:2211
void stick_leaves_into_vector(Vector< Tree * > &)
Traverse tree and stick pointers to leaf "nodes" (only) into Vector.
Definition: tree.cc:255
void add_external_halo_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add external halo element whose non-halo counterpart is held on processor p to this Mesh...
Definition: mesh.h:1949
SolidNode * boundary_node_pt(const unsigned &b, const unsigned &n)
Return n-th SolidNodes on b-th boundary.
Definition: mesh.h:2313
std::set< int > external_halo_proc()
Return the set of processors that hold external halo nodes. This is required to avoid having to pass ...
Definition: mesh.h:2160
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
unsigned ndof_types() const
Return number of dof types in mesh.
Definition: mesh.cc:8415
virtual void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Plot the error when compared against a given exact solution . Also calculates the norm of the error a...
Definition: elements.h:2992
void delete_all_external_storage()
Wipe the storage for all externally-based elements.
Definition: mesh.cc:8836
OomphCommunicator * Comm_pt
Pointer to communicator – set to NULL if mesh is not distributed.
Definition: mesh.h:130
void set_lagrangian_nodal_coordinates()
Make the current configuration the undeformed one by setting the nodal Lagrangian coordinates to thei...
Definition: mesh.cc:9196
virtual void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Plot error when compared against a given exact solution. Also returns the norm of the error and that ...
Definition: mesh.h:1096
bool is_halo() const
Is this element a halo?
Definition: elements.h:1141
A general mesh class.
Definition: mesh.h:74
void set_external_haloed_node_pt(const unsigned &p, const Vector< Node * > &external_haloed_node_pt)
Set vector of external haloed node in this Mesh whose halo external counterpart is held on processor ...
Definition: mesh.h:2151
GeneralisedElement *& external_haloed_element_pt(const unsigned &p, const unsigned &e)
Access fct to the e-th external haloed element in this Mesh whose non-halo counterpart is held on pro...
Definition: mesh.h:1985
IC problem for an elastic body discretised on a given (sub)-mesh. We switch the elements' residuals a...
void remove_boundary_nodes()
Clear all pointers to boundary nodes.
Definition: mesh.cc:208
Node * halo_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th halo node in this Mesh whose non-halo counterpart is held on processor p...
Definition: mesh.h:1616
SolidMesh(const Vector< SolidMesh * > &sub_mesh_pt)
Constructor builds combined mesh from the meshes specified. Note: This simply merges the meshes' elem...
Definition: mesh.h:2276
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types)...
virtual void create_shared_boundaries(OomphCommunicator *comm_pt, const Vector< unsigned > &element_domain, const Vector< GeneralisedElement * > &backed_up_el_pt, const Vector< FiniteElement * > &backed_up_f_el_pt, std::map< Data *, std::set< unsigned > > &processors_associated_with_data, const bool &overrule_keep_as_halo_element_status)
Creates the shared boundaries, only used in unstructured meshes In this case with the "TriangleMesh" ...
Definition: mesh.h:2173
void setup_shared_node_scheme()
Setup shared node scheme.
Definition: mesh.cc:2455
virtual void read(std::ifstream &restart_file)
Read solution from restart file.
Definition: mesh.cc:1087
void(FiniteElement::* UnsteadyExactSolutionFctPt)(const double &time, const Vector< double > &x, Vector< double > &soln)
Typedef for function pointer to function that computes unsteady exact solution.
Definition: mesh.h:241
void write_pvd_header(std::ofstream &pvd_file)
Write the pvd file header.
Definition: mesh.cc:9250
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:57