Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
oomph::TreeBasedRefineableMeshBase Class Referenceabstract

Base class for tree-based refineable meshes. More...

#include <refineable_mesh.h>

+ Inheritance diagram for oomph::TreeBasedRefineableMeshBase:

Public Member Functions

 TreeBasedRefineableMeshBase ()
 Constructor. More...
 
 TreeBasedRefineableMeshBase (const TreeBasedRefineableMeshBase &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const TreeBasedRefineableMeshBase &)=delete
 Broken assignment operator. More...
 
virtual ~TreeBasedRefineableMeshBase ()
 Empty Destructor: More...
 
void adapt (const Vector< double > &elemental_error)
 
void p_adapt (const Vector< double > &elemental_error)
 
void refine_uniformly (DocInfo &doc_info)
 Refine mesh uniformly and doc process. More...
 
void refine_uniformly ()
 Refine mesh uniformly. More...
 
void p_refine_uniformly (DocInfo &doc_info)
 p-refine mesh uniformly and doc process More...
 
void p_refine_uniformly ()
 p-refine mesh uniformly More...
 
unsigned unrefine_uniformly ()
 
void p_unrefine_uniformly (DocInfo &doc_info)
 p-unrefine mesh uniformly More...
 
virtual void setup_tree_forest ()=0
 Set up the tree forest associated with the Mesh (if any) More...
 
TreeForest * forest_pt ()
 Return pointer to the Forest represenation of the mesh. More...
 
void doc_adaptivity_targets (std::ostream &outfile)
 Doc the targets for mesh adaptation. More...
 
unsigned & max_refinement_level ()
 Access fct for max. permissible refinement level (relative to base mesh) More...
 
unsigned & min_refinement_level ()
 Access fct for min. permissible refinement level (relative to base mesh) More...
 
unsigned & max_p_refinement_level ()
 
unsigned & min_p_refinement_level ()
 
virtual void adapt_mesh (DocInfo &doc_info)
 
virtual void adapt_mesh ()
 
void p_adapt_mesh (DocInfo &doc_info)
 
void p_adapt_mesh ()
 
virtual void refine_selected_elements (const Vector< unsigned > &elements_to_be_refined)
 
virtual void refine_selected_elements (const Vector< RefineableElement * > &elements_to_be_refined)
 
void p_refine_selected_elements (const Vector< unsigned > &elements_to_be_refined)
 p-refine mesh by refining the elements identified by their numbers. More...
 
void p_refine_selected_elements (const Vector< PRefineableElement * > &elements_to_be_refined_pt)
 p-refine mesh by refining the elements identified by their pointers. More...
 
virtual void refine_base_mesh_as_in_reference_mesh (TreeBasedRefineableMeshBase *const &ref_mesh_pt)
 Refine to same degree as the reference mesh. More...
 
virtual bool refine_base_mesh_as_in_reference_mesh_minus_one (TreeBasedRefineableMeshBase *const &ref_mesh_pt)
 
virtual void refine_as_in_reference_mesh (TreeBasedRefineableMeshBase *const &ref_mesh_pt)
 
virtual void get_refinement_levels (unsigned &min_refinement_level, unsigned &max_refinement_level)
 Get max/min refinement levels in mesh. More...
 
virtual void get_elements_at_refinement_level (unsigned &refinement_level, Vector< RefineableElement * > &level_elements)
 
virtual void get_refinement_pattern (Vector< Vector< unsigned >> &to_be_refined)
 
void refine_base_mesh (Vector< Vector< unsigned >> &to_be_refined)
 Refine base mesh according to specified refinement pattern. More...
 
virtual void refine (std::ifstream &restart_file)
 Refine mesh according to refinement pattern in restart file. More...
 
virtual void dump_refinement (std::ostream &outfile)
 Dump refinement pattern to allow for rebuild. More...
 
virtual void read_refinement (std::ifstream &restart_file, Vector< Vector< unsigned >> &to_be_refined)
 Read refinement pattern to allow for rebuild. More...
 
unsigned uniform_refinement_level_when_pruned () const
 
unsigned & uniform_refinement_level_when_pruned ()
 Level to which the mesh was uniformly refined when it was pruned. More...
 
- Public Member Functions inherited from oomph::RefineableMeshBase
bool adapt_flag ()
 
 RefineableMeshBase ()
 
 RefineableMeshBase (const RefineableMeshBase &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const RefineableMeshBase &)=delete
 Broken assignment operator. More...
 
virtual ~RefineableMeshBase ()
 Empty Destructor: More...
 
unsigned nrefined ()
 Access fct for number of elements that were refined. More...
 
unsigned nunrefined ()
 Access fct for number of elements that were unrefined. More...
 
unsigned & nrefinement_overruled ()
 
unsigned & max_keep_unrefined ()
 
ErrorEstimator *& spatial_error_estimator_pt ()
 Access to spatial error estimator. More...
 
ErrorEstimator * spatial_error_estimator_pt () const
 Access to spatial error estimator (const version. More...
 
double & min_permitted_error ()
 
double & max_permitted_error ()
 
double & min_error ()
 
double & max_error ()
 
DocInfo *& doc_info_pt ()
 Access fct for pointer to DocInfo. More...
 
void enable_adaptation ()
 Enable adaptation. More...
 
void disable_adaptation ()
 Disable adaptation. More...
 
void enable_p_adaptation ()
 Enable adaptation. More...
 
void disable_p_adaptation ()
 Disable adaptation. More...
 
void enable_additional_synchronisation_of_hanging_nodes ()
 Enable additional synchronisation of hanging nodes. More...
 
void disable_additional_synchronisation_of_hanging_nodes ()
 Disable additional synchronisation of hanging nodes. More...
 
bool is_adaptation_enabled () const
 Return whether the mesh is to be adapted. More...
 
bool is_p_adaptation_enabled () const
 Return whether the mesh is to be adapted. More...
 
bool is_additional_synchronisation_of_hanging_nodes_disabled () const
 Return whether additional synchronisation is enabled. More...
 
DocInfo doc_info ()
 Access fct for DocInfo. More...
 
void p_unrefine_uniformly (DocInfo &doc_info)
 p-unrefine mesh uniformly More...
 
- Public Member Functions inherited from oomph::Mesh
 Mesh ()
 Default constructor. More...
 
 Mesh (const Vector< Mesh * > &sub_mesh_pt)
 
void merge_meshes (const Vector< Mesh * > &sub_mesh_pt)
 
virtual void setup_boundary_element_info ()
 
virtual void setup_boundary_element_info (std::ostream &outfile)
 
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. More...
 
template<class BULK_ELEMENT >
void doc_boundary_coordinates (const unsigned &b, std::ofstream &the_file)
 
virtual void scale_mesh (const double &factor)
 
 Mesh (const Mesh &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const Mesh &)=delete
 Broken assignment operator. More...
 
virtual ~Mesh ()
 Virtual Destructor to clean up all memory. More...
 
void flush_element_and_node_storage ()
 
void flush_element_storage ()
 
void flush_node_storage ()
 
Node *& node_pt (const unsigned long &n)
 Return pointer to global node n. More...
 
Node * node_pt (const unsigned long &n) const
 Return pointer to global node n (const version) More...
 
GeneralisedElement *& element_pt (const unsigned long &e)
 Return pointer to element e. More...
 
GeneralisedElement * element_pt (const unsigned long &e) const
 Return pointer to element e (const version) More...
 
const Vector< GeneralisedElement * > & element_pt () const
 Return reference to the Vector of elements. More...
 
Vector< GeneralisedElement * > & element_pt ()
 Return reference to the Vector of elements. More...
 
FiniteElement * finite_element_pt (const unsigned &e) const
 
Node *& boundary_node_pt (const unsigned &b, const unsigned &n)
 Return pointer to node n on boundary b. More...
 
Node * boundary_node_pt (const unsigned &b, const unsigned &n) const
 Return pointer to node n on boundary b. More...
 
void set_nboundary (const unsigned &nbound)
 Set the number of boundaries in the mesh. More...
 
void remove_boundary_nodes ()
 Clear all pointers to boundary nodes. More...
 
void remove_boundary_nodes (const unsigned &b)
 
void remove_boundary_node (const unsigned &b, Node *const &node_pt)
 Remove a node from the boundary b. More...
 
void add_boundary_node (const unsigned &b, Node *const &node_pt)
 Add a (pointer to) a node to the b-th boundary. More...
 
void copy_boundary_node_data_from_nodes ()
 
bool boundary_coordinate_exists (const unsigned &i) const
 Indicate whether the i-th boundary has an intrinsic coordinate. More...
 
unsigned long nelement () const
 Return number of elements in the mesh. More...
 
unsigned long nnode () const
 Return number of nodes in the mesh. More...
 
unsigned ndof_types () const
 Return number of dof types in mesh. More...
 
unsigned elemental_dimension () const
 Return number of elemental dimension in mesh. More...
 
unsigned nodal_dimension () const
 Return number of nodal dimension in mesh. More...
 
void add_node_pt (Node *const &node_pt)
 Add a (pointer to a) node to the mesh. More...
 
void add_element_pt (GeneralisedElement *const &element_pt)
 Add a (pointer to) an element to the mesh. More...
 
virtual void node_update (const bool &update_all_solid_nodes=false)
 
virtual void reorder_nodes (const bool &use_old_ordering=true)
 
virtual void get_node_reordering (Vector< Node * > &reordering, const bool &use_old_ordering=true) const
 
template<class BULK_ELEMENT , template< class > class FACE_ELEMENT>
void build_face_mesh (const unsigned &b, Mesh *const &face_mesh_pt)
 
unsigned self_test ()
 Self-test: Check elements and nodes. Return 0 for OK. More...
 
void max_and_min_element_size (double &max_size, double &min_size)
 
double total_size ()
 
void check_inverted_elements (bool &mesh_has_inverted_elements, std::ofstream &inverted_element_file)
 
void check_inverted_elements (bool &mesh_has_inverted_elements)
 
unsigned check_for_repeated_nodes (const double &epsilon=1.0e-12)
 
Vector< Node * > prune_dead_nodes ()
 
unsigned nboundary () const
 Return number of boundaries. More...
 
unsigned long nboundary_node (const unsigned &ibound) const
 Return number of nodes on a particular boundary. More...
 
FiniteElement * boundary_element_pt (const unsigned &b, const unsigned &e) const
 Return pointer to e-th finite element on boundary b. More...
 
Node * get_some_non_boundary_node () const
 
unsigned nboundary_element (const unsigned &b) const
 Return number of finite elements that are adjacent to boundary b. More...
 
int face_index_at_boundary (const unsigned &b, const unsigned &e) const
 
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. More...
 
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. More...
 
virtual void read (std::ifstream &restart_file)
 Read solution from restart file. More...
 
void output_paraview (std::ofstream &file_out, const unsigned &nplot) const
 
void output_fct_paraview (std::ofstream &file_out, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
 
void output_fct_paraview (std::ofstream &file_out, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
void output (std::ostream &outfile)
 Output for all elements. More...
 
void output (std::ostream &outfile, const unsigned &n_plot)
 Output at f(n_plot) points in each element. More...
 
void output (FILE *file_pt)
 Output for all elements (C-style output) More...
 
void output (FILE *file_pt, const unsigned &nplot)
 Output at f(n_plot) points in each element (C-style output) More...
 
void output (const std::string &output_filename)
 Output for all elements. More...
 
void output (const std::string &output_filename, const unsigned &n_plot)
 Output at f(n_plot) points in each element. More...
 
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. More...
 
void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt)
 
void output_boundaries (std::ostream &outfile)
 Output the nodes on the boundaries (into separate tecplot zones) More...
 
void output_boundaries (const std::string &output_filename)
 
void assign_initial_values_impulsive ()
 Assign initial values for an impulsive start. More...
 
void shift_time_values ()
 
void calculate_predictions ()
 
void set_nodal_and_elemental_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
virtual void set_mesh_level_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void set_consistent_pinned_values_for_continuation (ContinuationStorageScheme *const &continuation_stepper_pt)
 Set consistent values for pinned data in continuation. More...
 
bool does_pointer_correspond_to_mesh_data (double *const &parameter_pt)
 Does the double pointer correspond to any mesh data. More...
 
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. More...
 
void set_elemental_internal_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
virtual void compute_norm (double &norm)
 
virtual void compute_norm (Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 Returns the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
bool is_mesh_distributed () const
 Boolean to indicate if Mesh has been distributed. More...
 
OomphCommunicator * communicator_pt () const
 
void delete_all_external_storage ()
 Wipe the storage for all externally-based elements. More...
 

Protected Member Functions

virtual void split_elements_if_required ()=0
 
virtual void p_refine_elements_if_required ()=0
 
void complete_hanging_nodes (const int &ncont_interpolated_values)
 Complete the hanging node scheme recursively. More...
 
void complete_hanging_nodes_recursively (Node *&nod_pt, Vector< Node * > &master_nodes, Vector< double > &hang_weights, const int &ival)
 Auxiliary routine for recursive hanging node completion. More...
 
- Protected Member Functions inherited from oomph::Mesh
unsigned long assign_global_eqn_numbers (Vector< double * > &Dof_pt)
 Assign (global) equation numbers to the nodes. More...
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 Assign local equation numbers in all elements. More...
 
void convert_to_boundary_node (Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
 
void convert_to_boundary_node (Node *&node_pt)
 

Protected Attributes

unsigned Uniform_refinement_level_when_pruned
 Level to which the mesh was uniformly refined when it was pruned. More...
 
unsigned Max_refinement_level
 Max. permissible refinement level (relative to base mesh) More...
 
unsigned Min_refinement_level
 Min. permissible refinement level (relative to base mesh) More...
 
unsigned Max_p_refinement_level
 Max. permissible p-refinement level (relative to base mesh) More...
 
unsigned Min_p_refinement_level
 Min. permissible p-refinement level (relative to base mesh) More...
 
TreeForest * Forest_pt
 Forest representation of the mesh. More...
 
- Protected Attributes inherited from oomph::RefineableMeshBase
ErrorEstimator * Spatial_error_estimator_pt
 Pointer to spatial error estimator. More...
 
double Max_permitted_error
 Max. error (i.e. split elements if their error is larger) More...
 
double Min_permitted_error
 Min. error (i.e. (try to) merge elements if their error is smaller) More...
 
double Min_error
 Min.actual error. More...
 
double Max_error
 Max. actual error. More...
 
unsigned Nrefined
 Stats: Number of elements that were refined. More...
 
unsigned Nunrefined
 Stats: Number of elements that were unrefined. More...
 
bool Adapt_flag
 Flag that requests adaptation. More...
 
bool P_adapt_flag
 Flag that requests p-adaptation. More...
 
bool Additional_synchronisation_of_hanging_nodes_not_required
 Flag that disables additional synchronisation of hanging nodes. More...
 
DocInfo * Doc_info_pt
 Pointer to DocInfo. More...
 
unsigned Max_keep_unrefined
 
unsigned Nrefinement_overruled
 
- Protected Attributes inherited from oomph::Mesh
Vector< Vector< Node * > > Boundary_node_pt
 
bool Lookup_for_elements_next_boundary_is_setup
 
Vector< Vector< FiniteElement * > > Boundary_element_pt
 
Vector< Vector< int > > Face_index_at_boundary
 
Vector< Node * > Node_pt
 Vector of pointers to nodes. More...
 
Vector< GeneralisedElement * > Element_pt
 Vector of pointers to generalised elements. More...
 
std::vector< bool > Boundary_coordinate_exists
 

Additional Inherited Members

- Public Types inherited from oomph::Mesh
typedef void(FiniteElement::* SteadyExactSolutionFctPt) (const Vector< double > &x, Vector< double > &soln)
 
typedef void(FiniteElement::* UnsteadyExactSolutionFctPt) (const double &time, const Vector< double > &x, Vector< double > &soln)
 
- Static Public Attributes inherited from oomph::Mesh
static Steady< 0 > Default_TimeStepper
 The Steady Timestepper. More...
 
static bool Suppress_warning_about_empty_mesh_level_time_stepper_function
 Static boolean flag to control warning about mesh level timesteppers. More...
 

Detailed Description

Base class for tree-based refineable meshes.

///////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////

Constructor & Destructor Documentation

◆ TreeBasedRefineableMeshBase() [1/2]

oomph::TreeBasedRefineableMeshBase::TreeBasedRefineableMeshBase ( )
inline

Constructor.

380  {
381  // Max/min refinement levels
382  Max_refinement_level = 5;
383  Min_refinement_level = 0;
384 
385  // Max/min p-refinement levels
386  Max_p_refinement_level = 7;
387  Min_p_refinement_level = 2;
388 
389  // Stats
390  Nrefined = 0;
391  Nunrefined = 0;
392 
393  // Where do I write the documatation of the refinement process?
394  Doc_info_pt = 0;
395 
396  // Initialise the forest pointer to NULL
397  Forest_pt = 0;
398 
399  // Mesh hasn't been pruned yet
400  Uniform_refinement_level_when_pruned = 0;
401  }
oomph::RefineableMeshBase::Nrefined
unsigned Nrefined
Stats: Number of elements that were refined.
Definition: refineable_mesh.h:338
oomph::RefineableMeshBase::Doc_info_pt
DocInfo * Doc_info_pt
Pointer to DocInfo.
Definition: refineable_mesh.h:353
oomph::RefineableMeshBase::Nunrefined
unsigned Nunrefined
Stats: Number of elements that were unrefined.
Definition: refineable_mesh.h:341
oomph::TreeBasedRefineableMeshBase::Uniform_refinement_level_when_pruned
unsigned Uniform_refinement_level_when_pruned
Level to which the mesh was uniformly refined when it was pruned.
Definition: refineable_mesh.h:753
oomph::TreeBasedRefineableMeshBase::Max_p_refinement_level
unsigned Max_p_refinement_level
Max. permissible p-refinement level (relative to base mesh)
Definition: refineable_mesh.h:762
oomph::TreeBasedRefineableMeshBase::Forest_pt
TreeForest * Forest_pt
Forest representation of the mesh.
Definition: refineable_mesh.h:768
oomph::TreeBasedRefineableMeshBase::Min_p_refinement_level
unsigned Min_p_refinement_level
Min. permissible p-refinement level (relative to base mesh)
Definition: refineable_mesh.h:765
oomph::TreeBasedRefineableMeshBase::Max_refinement_level
unsigned Max_refinement_level
Max. permissible refinement level (relative to base mesh)
Definition: refineable_mesh.h:756
oomph::TreeBasedRefineableMeshBase::Min_refinement_level
unsigned Min_refinement_level
Min. permissible refinement level (relative to base mesh)
Definition: refineable_mesh.h:759

References oomph::RefineableMeshBase::Doc_info_pt, Forest_pt, Max_p_refinement_level, Max_refinement_level, Min_p_refinement_level, Min_refinement_level, oomph::RefineableMeshBase::Nrefined, oomph::RefineableMeshBase::Nunrefined, and Uniform_refinement_level_when_pruned.

◆ TreeBasedRefineableMeshBase() [2/2]

oomph::TreeBasedRefineableMeshBase::TreeBasedRefineableMeshBase ( const TreeBasedRefineableMeshBase &  dummy)
delete

Broken copy constructor.

◆ ~TreeBasedRefineableMeshBase()

virtual oomph::TreeBasedRefineableMeshBase::~TreeBasedRefineableMeshBase ( )
inlinevirtual

Empty Destructor:

413  {
414  // Kill the forest if there is one
415  if (Forest_pt != 0)
416  {
417  delete Forest_pt;
418  Forest_pt = 0;
419  }
420  }

References Forest_pt.

Member Function Documentation

◆ adapt()

void oomph::TreeBasedRefineableMeshBase::adapt ( const Vector< double > &  elemental_error)
virtual

Adapt mesh: Refine elements whose error is lager than err_max and (try to) unrefine those whose error is smaller than err_min

Do adaptive refinement for mesh.

  • Pass Vector of error estimates for all elements.
  • Refine those whose errors exceeds the threshold
  • (Try to) unrefine those whose errors is less than threshold (only possible if the three brothers also want to be unrefined, of course.)
  • Update the nodal positions in the whole lot
  • Store # of refined/unrefined elements.
  • Doc refinement process (if required)

Implements oomph::RefineableMeshBase.

Reimplemented in oomph::RefineableQuadFromTriangleMesh< ELEMENT >.

308  {
309  // Set the refinement tolerance to be the max permissible error
310  double refine_tol = max_permitted_error();
311 
312  // Set the unrefinement tolerance to be the min permissible error
313  double unrefine_tol = min_permitted_error();
314 
315  // Setup doc info
316  DocInfo local_doc_info;
317  if (doc_info_pt() == 0)
318  {
319  local_doc_info.disable_doc();
320  }
321  else
322  {
323  local_doc_info = doc_info();
324  }
325 
326 
327  // Check that the errors make sense
328  if (refine_tol <= unrefine_tol)
329  {
330  std::ostringstream error_stream;
331  error_stream << "Refinement tolerance <= Unrefinement tolerance"
332  << refine_tol << " " << unrefine_tol << std::endl
333  << "doesn't make sense and will almost certainly crash"
334  << std::endl
335  << "this beautiful code!" << std::endl;
336 
337  throw OomphLibError(
338  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
339  }
340 
341 
342  // Select elements for refinement and unrefinement
343  //================================================
344  // Reset counter for number of elements that would like to be
345  // refined further but can't
346  nrefinement_overruled() = 0;
347 
348  // Note: Yes, this needs to be a map because we'll have to check
349  // the refinement wishes of brothers (who we only access via pointers)
350  std::map<RefineableElement*, bool> wants_to_be_unrefined;
351 
352  // Initialise a variable to store the number of elements for refinment
353  unsigned n_refine = 0;
354 
355  // Loop over all elements and mark them according to the error criterion
356  unsigned long Nelement = this->nelement();
357  for (unsigned long e = 0; e < Nelement; e++)
358  {
359  // (Cast) pointer to the element
360  RefineableElement* el_pt =
361  dynamic_cast<RefineableElement*>(this->element_pt(e));
362 
363  // Initially element is not to be refined
364  el_pt->deselect_for_refinement();
365 
366  // If the element error exceeds the threshold ...
367  if (elemental_error[e] > refine_tol)
368  {
369  // ... and its refinement level is less than the maximum desired level
370  // mark is to be refined
371  if ((el_pt->refinement_is_enabled()) &&
372  (el_pt->refinement_level() < max_refinement_level()))
373  {
374  el_pt->select_for_refinement();
375  n_refine++;
376  }
377  // ... otherwise mark it as having been over-ruled
378  else
379  {
380  nrefinement_overruled() += 1;
381  }
382  }
383 
384  // Now worry about unrefinement (first pass):
385 
386  // Is my error too small AND do I have a father?
387  if ((elemental_error[e] < unrefine_tol) &&
388  (el_pt->tree_pt()->father_pt() != 0))
389  {
390  // Flag to indicate whether to unrefine
391  bool unrefine = true;
392  unsigned n_sons = el_pt->tree_pt()->father_pt()->nsons();
393 
394  // Are all brothers leaf nodes?
395  for (unsigned ison = 0; ison < n_sons; ison++)
396  {
397  // (At least) one brother is not a leaf: end of story; we're not doing
398  // it (= the unrefinement)
399  if (!(el_pt->tree_pt()->father_pt()->son_pt(ison)->is_leaf()))
400  {
401  unrefine = false;
402  }
403  }
404 
405  // Don't allow unrefinement of elements that would become larger
406  // than the minimum legal refinement level
407  if (el_pt->refinement_level() - 1 < min_refinement_level())
408  {
409  unrefine = false;
410  }
411 
412  // So, all things considered, is the element eligbible for refinement?
413  if (unrefine)
414  {
415  wants_to_be_unrefined[el_pt] = true;
416  }
417  else
418  {
419  wants_to_be_unrefined[el_pt] = false;
420  }
421  }
422  }
423 
424  oomph_info << " \n Number of elements to be refined: " << n_refine
425  << std::endl;
426  oomph_info << " \n Number of elements whose refinement was overruled: "
427  << nrefinement_overruled() << std::endl;
428 
429  // Second pass for unrefinement --- an element cannot be unrefined unless
430  // all brothers want to be unrefined.
431  // Loop over all elements again and let the first set of sons check if their
432  // brothers also want to be unrefined
433  unsigned n_unrefine = 0;
434  for (unsigned long e = 0; e < Nelement; e++)
435  {
436  //(Cast) pointer to the element
437  RefineableElement* el_pt =
438  dynamic_cast<RefineableElement*>(this->element_pt(e));
439 
440  // hierher: This is a bit naughty... We want to put the
441  // first son in charge -- the statement below assumes (correctly) that the
442  // enumeration of all (!) trees starts with son types.
443  // This is correct for oc and quadtrees but will bite us if we
444  // ever introduce other trees if/when we accidentally break this
445  // tacit assumption. Not sure what to do about it for
446  // now other than leaving it hierher-ed...
447  if (el_pt->tree_pt()->son_type() == OcTreeNames::LDB)
448  {
449  // Do all sons want to be unrefined?
450  bool unrefine = true;
451  unsigned n_sons = el_pt->tree_pt()->father_pt()->nsons();
452  for (unsigned ison = 0; ison < n_sons; ison++)
453  {
454  if (!(wants_to_be_unrefined[dynamic_cast<RefineableElement*>(
455  el_pt->tree_pt()->father_pt()->son_pt(ison)->object_pt())]))
456  {
457  // One guy isn't cooperating and spoils the party.
458  unrefine = false;
459  }
460  }
461 
462  // Tell father that his sons need to be merged
463  if (unrefine)
464  {
465  el_pt->tree_pt()
466  ->father_pt()
467  ->object_pt()
468  ->select_sons_for_unrefinement();
469  n_unrefine += n_sons;
470  }
471  // Otherwise mark the sons as not to be touched
472  else
473  {
474  el_pt->tree_pt()
475  ->father_pt()
476  ->object_pt()
477  ->deselect_sons_for_unrefinement();
478  }
479  }
480  }
481  oomph_info << " \n Number of elements to be merged : " << n_unrefine
482  << std::endl
483  << std::endl;
484 
485 
486  // Now do the actual mesh adaptation
487  //---------------------------------
488 
489  // Check whether its worth our while
490  // Either some elements want to be refined,
491  // or the number that want to be unrefined are greater than the
492  // specified tolerance
493 
494  // In a parallel job, it is possible that one process may not have
495  // any elements to refine, BUT a neighbouring process may refine an
496  // element which changes the hanging status of a node that is on
497  // both processes (i.e. a halo(ed) node). To get around this issue,
498  // ALL processes need to call adapt_mesh if ANY refinement is to
499  // take place anywhere.
500 
501  unsigned total_n_refine = 0;
502 #ifdef OOMPH_HAS_MPI
503  // Sum n_refine across all processors
504  if (this->is_mesh_distributed())
505  {
506  MPI_Allreduce(&n_refine,
507  &total_n_refine,
508  1,
509  MPI_UNSIGNED,
510  MPI_SUM,
511  Comm_pt->mpi_comm());
512  }
513  else
514  {
515  total_n_refine = n_refine;
516  }
517 #else
518  total_n_refine = n_refine;
519 #endif
520 
521  // There may be some issues with unrefinement too, but I have not
522  // been able to come up with an example (either in my head or in a
523  // particular problem) where anything has arisen. I can see that
524  // there may be an issue if n_unrefine differs across processes so
525  // that (total_n_unrefine > max_keep_unrefined()) on some but not
526  // all processes. I haven't seen any examples of this yet so the
527  // following code may or may not work! (Andy, 06/03/08)
528 
529  unsigned total_n_unrefine = 0;
530 #ifdef OOMPH_HAS_MPI
531  // Sum n_unrefine across all processors
532  if (this->is_mesh_distributed())
533  {
534  MPI_Allreduce(&n_unrefine,
535  &total_n_unrefine,
536  1,
537  MPI_UNSIGNED,
538  MPI_SUM,
539  Comm_pt->mpi_comm());
540  }
541  else
542  {
543  total_n_unrefine = n_unrefine;
544  }
545 #else
546  total_n_unrefine = n_unrefine;
547 #endif
548 
549  oomph_info << "---> " << total_n_refine << " elements to be refined, and "
550  << total_n_unrefine << " to be unrefined, in total.\n"
551  << std::endl;
552 
553  if ((total_n_refine > 0) || (total_n_unrefine > max_keep_unrefined()))
554  {
555 #ifdef PARANOID
556 #ifdef OOMPH_HAS_MPI
557 
558 
559  // Sanity check: Each processor checks if the enforced unrefinement of
560  // its haloed element is matched by enforced unrefinement of the
561  // corresponding halo elements on the other processors.
562  if (this->is_mesh_distributed())
563  {
564  // Store number of processors and current process
565  MPI_Status status;
566  int n_proc = Comm_pt->nproc();
567  int my_rank = Comm_pt->my_rank();
568 
569  // Loop over processes: Each processor sends unrefinement pattern
570  // for halo elements with processor d to processor d where it's
571  // compared against the unrefinement pattern for the corresponding
572  // haloed elements
573  for (int d = 0; d < n_proc; d++)
574  {
575  // No halo with self: Send unrefinement info to proc d
576  if (d != my_rank)
577  {
578  // Get the vector of halo elements whose non-halo counterpart
579  // are on processor d
580  Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(d));
581 
582  // Create vector containing (0)1 to indicate that
583  // halo element is (not) to be unrefined
584  unsigned nhalo = halo_elem_pt.size();
585  Vector<int> halo_to_be_unrefined(nhalo, 0);
586  for (unsigned e = 0; e < nhalo; e++)
587  {
588  if (dynamic_cast<RefineableElement*>(halo_elem_pt[e])
589  ->sons_to_be_unrefined())
590  {
591  halo_to_be_unrefined[e] = 1;
592  }
593  }
594 
595  // Trap the case when there are no halo elements
596  // so that we don't get a segfault in the MPI send
597  if (nhalo > 0)
598  {
599  // Send it across to proc d
600  MPI_Send(&halo_to_be_unrefined[0],
601  nhalo,
602  MPI_INT,
603  d,
604  0,
605  Comm_pt->mpi_comm());
606  }
607  }
608  // else (d=my_rank): Receive unrefinement pattern from all
609  // other processors (dd)
610  else
611  {
612  // Loop over other processors
613  for (int dd = 0; dd < n_proc; dd++)
614  {
615  // No halo with yourself
616  if (dd != d)
617  {
618  // Get the vector of haloed elements on current processor
619  // with processor dd
620  Vector<GeneralisedElement*> haloed_elem_pt(
621  this->haloed_element_pt(dd));
622 
623  // Ask processor dd to send vector containing (0)1 for
624  // halo element with current processor to be (not)unrefined
625  unsigned nhaloed = haloed_elem_pt.size();
626  Vector<int> halo_to_be_unrefined(nhaloed);
627  // Trap to catch the case that there are no haloed elements
628  if (nhaloed > 0)
629  {
630  // Receive unrefinement pattern of haloes from proc dd
631  MPI_Recv(&halo_to_be_unrefined[0],
632  nhaloed,
633  MPI_INT,
634  dd,
635  0,
636  Comm_pt->mpi_comm(),
637  &status);
638  }
639 
640  // Check it
641  for (unsigned e = 0; e < nhaloed; e++)
642  {
643  if (((halo_to_be_unrefined[e] == 0) &&
644  (dynamic_cast<RefineableElement*>(haloed_elem_pt[e])
645  ->sons_to_be_unrefined())) ||
646  ((halo_to_be_unrefined[e] == 1) &&
647  (!dynamic_cast<RefineableElement*>(haloed_elem_pt[e])
648  ->sons_to_be_unrefined())))
649  {
650  std::ostringstream error_message;
651  error_message
652  << "Error in unrefinement: \n"
653  << "Haloed element: " << e << " on proc " << my_rank
654  << " \n"
655  << "wants to be unrefined whereas its halo counterpart "
656  "on\n"
657  << "proc " << dd << " doesn't (or vice versa)...\n"
658  << "This is most likely because the error estimator\n"
659  << "has not assigned the same errors to halo and haloed\n"
660  << "elements -- it ought to!\n";
661  throw OomphLibError(error_message.str(),
662  OOMPH_CURRENT_FUNCTION,
663  OOMPH_EXCEPTION_LOCATION);
664  }
665  }
666  }
667  }
668  }
669  }
670 
671 
672  // Loop over processes: Each processor sends refinement pattern
673  // for halo elements with processor d to processor d where it's
674  // compared against the refinement pattern for the corresponding
675  // haloed elements
676  for (int d = 0; d < n_proc; d++)
677  {
678  // No halo with self: Send refinement info to proc d
679  if (d != my_rank)
680  {
681  // Get the vector of halo elements whose non-halo counterpart
682  // are on processor d
683  Vector<GeneralisedElement*> halo_elem_pt(this->halo_element_pt(d));
684 
685  // Create vector containing (0)1 to indicate that
686  // halo element is (not) to be refined
687  unsigned nhalo = halo_elem_pt.size();
688  Vector<int> halo_to_be_refined(nhalo, 0);
689  for (unsigned e = 0; e < nhalo; e++)
690  {
691  if (dynamic_cast<RefineableElement*>(halo_elem_pt[e])
692  ->to_be_refined())
693  {
694  halo_to_be_refined[e] = 1;
695  }
696  }
697 
698  // Trap the case when there are no halo elements
699  // so that we don't get a segfault in the MPI send
700  if (nhalo > 0)
701  {
702  // Send it across to proc d
703  MPI_Send(&halo_to_be_refined[0],
704  nhalo,
705  MPI_INT,
706  d,
707  0,
708  Comm_pt->mpi_comm());
709  }
710  }
711  // else (d=my_rank): Receive refinement pattern from all
712  // other processors (dd)
713  else
714  {
715  // Loop over other processors
716  for (int dd = 0; dd < n_proc; dd++)
717  {
718  // No halo with yourself
719  if (dd != d)
720  {
721  // Get the vector of haloed elements on current processor
722  // with processor dd
723  Vector<GeneralisedElement*> haloed_elem_pt(
724  this->haloed_element_pt(dd));
725 
726  // Ask processor dd to send vector containing (0)1 for
727  // halo element with current processor to be (not)refined
728  unsigned nhaloed = haloed_elem_pt.size();
729  Vector<int> halo_to_be_refined(nhaloed);
730  // Trap to catch the case that there are no haloed elements
731  if (nhaloed > 0)
732  {
733  // Receive unrefinement pattern of haloes from proc dd
734  MPI_Recv(&halo_to_be_refined[0],
735  nhaloed,
736  MPI_INT,
737  dd,
738  0,
739  Comm_pt->mpi_comm(),
740  &status);
741  }
742 
743  // Check it
744  for (unsigned e = 0; e < nhaloed; e++)
745  {
746  if (((halo_to_be_refined[e] == 0) &&
747  (dynamic_cast<RefineableElement*>(haloed_elem_pt[e])
748  ->to_be_refined())) ||
749  ((halo_to_be_refined[e] == 1) &&
750  (!dynamic_cast<RefineableElement*>(haloed_elem_pt[e])
751  ->to_be_refined())))
752  {
753  std::ostringstream error_message;
754  error_message
755  << "Error in refinement: \n"
756  << "Haloed element: " << e << " on proc " << my_rank
757  << " \n"
758  << "wants to be refined whereas its halo counterpart on\n"
759  << "proc " << dd << " doesn't (or vice versa)...\n"
760  << "This is most likely because the error estimator\n"
761  << "has not assigned the same errors to halo and haloed\n"
762  << "elements -- it ought to!\n";
763  throw OomphLibError(error_message.str(),
764  OOMPH_CURRENT_FUNCTION,
765  OOMPH_EXCEPTION_LOCATION);
766  }
767  }
768  }
769  }
770  }
771  }
772  }
773 
774 #endif
775 #endif
776 
777  // Perform the actual adaptation
778  adapt_mesh(local_doc_info);
779 
780  // The number of refineable elements is still local to each process
781  Nunrefined = n_unrefine;
782  Nrefined = n_refine;
783  }
784  // If not worthwhile, say so but still reorder nodes and kill
785  // external storage for consistency in parallel computations
786  else
787  {
788 #ifdef OOMPH_HAS_MPI
789  // Delete any external element storage - any interaction will still
790  // be set up on the fly again, so we need to get rid of old information.
791  // This particularly causes problems in multi-domain examples where
792  // we decide not to refine one of the meshes
793  this->delete_all_external_storage();
794 #endif
795 
796  // Reorder the nodes within the mesh's node vector
797  // to establish a standard ordering regardless of the sequence
798  // of mesh refinements -- this is required to allow dump/restart
799  // on refined meshes
800  this->reorder_nodes();
801 
802 #ifdef OOMPH_HAS_MPI
803 
804  // Now (re-)classify halo and haloed nodes and synchronise hanging
805  // nodes
806  // This is required in cases where delete_all_external_storage()
807  // made dependent nodes to external halo nodes nonhanging.
808  if (this->is_mesh_distributed())
809  {
810  DocInfo doc_info;
811  doc_info.disable_doc();
812  classify_halo_and_haloed_nodes(doc_info, doc_info.is_doc_enabled());
813  }
814 
815 #endif
816 
817  if (n_refine == 0)
818  {
819  oomph_info << " Not enough benefit in adapting mesh." << std::endl
820  << std::endl;
821  }
822  Nunrefined = 0;
823  Nrefined = 0;
824  }
825  }
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
oomph::DocInfo::is_doc_enabled
bool is_doc_enabled() const
Are we documenting?
Definition: oomph_utilities.h:548
oomph::DocInfo::disable_doc
void disable_doc()
Disable documentation.
Definition: oomph_utilities.h:542
oomph::Mesh::is_mesh_distributed
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1588
oomph::Mesh::reorder_nodes
virtual void reorder_nodes(const bool &use_old_ordering=true)
Definition: mesh.cc:508
oomph::Mesh::delete_all_external_storage
void delete_all_external_storage()
Wipe the storage for all externally-based elements.
Definition: mesh.cc:9190
oomph::Mesh::element_pt
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
Definition: mesh.h:460
oomph::Mesh::nelement
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
oomph::RefineableMeshBase::max_keep_unrefined
unsigned & max_keep_unrefined()
Definition: refineable_mesh.h:123
oomph::RefineableMeshBase::doc_info
DocInfo doc_info()
Access fct for DocInfo.
Definition: refineable_mesh.h:243
oomph::RefineableMeshBase::doc_info_pt
DocInfo *& doc_info_pt()
Access fct for pointer to DocInfo.
Definition: refineable_mesh.h:183
oomph::RefineableMeshBase::min_permitted_error
double & min_permitted_error()
Definition: refineable_mesh.h:156
oomph::RefineableMeshBase::max_permitted_error
double & max_permitted_error()
Definition: refineable_mesh.h:163
oomph::RefineableMeshBase::nrefinement_overruled
unsigned & nrefinement_overruled()
Definition: refineable_mesh.h:114
oomph::TreeBasedRefineableMeshBase::max_refinement_level
unsigned & max_refinement_level()
Access fct for max. permissible refinement level (relative to base mesh)
Definition: refineable_mesh.h:486
oomph::TreeBasedRefineableMeshBase::adapt_mesh
virtual void adapt_mesh()
Definition: refineable_mesh.h:517
oomph::TreeBasedRefineableMeshBase::min_refinement_level
unsigned & min_refinement_level()
Access fct for min. permissible refinement level (relative to base mesh)
Definition: refineable_mesh.h:492
oomph::OcTreeNames::LDB
@ LDB
Definition: octree.h:49
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References adapt_mesh(), oomph::Mesh::delete_all_external_storage(), oomph::RefineableElement::deselect_for_refinement(), oomph::RefineableElement::deselect_sons_for_unrefinement(), oomph::DocInfo::disable_doc(), oomph::RefineableMeshBase::doc_info(), oomph::RefineableMeshBase::doc_info_pt(), e(), oomph::Mesh::element_pt(), oomph::Tree::father_pt(), oomph::DocInfo::is_doc_enabled(), oomph::Tree::is_leaf(), oomph::Mesh::is_mesh_distributed(), oomph::OcTreeNames::LDB, oomph::RefineableMeshBase::max_keep_unrefined(), oomph::RefineableMeshBase::max_permitted_error(), max_refinement_level(), oomph::RefineableMeshBase::min_permitted_error(), min_refinement_level(), oomph::Mesh::nelement(), oomph::RefineableMeshBase::Nrefined, oomph::RefineableMeshBase::nrefinement_overruled(), oomph::Tree::nsons(), oomph::RefineableMeshBase::Nunrefined, oomph::Tree::object_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::RefineableElement::refinement_is_enabled(), oomph::RefineableElement::refinement_level(), oomph::Mesh::reorder_nodes(), oomph::RefineableElement::select_for_refinement(), oomph::RefineableElement::select_sons_for_unrefinement(), oomph::Tree::son_pt(), oomph::Tree::son_type(), and oomph::RefineableElement::tree_pt().

Referenced by unrefine_uniformly().

◆ adapt_mesh() [1/2]

virtual void oomph::TreeBasedRefineableMeshBase::adapt_mesh ( )
inlinevirtual

Perform the actual tree-based mesh adaptation. A simple wrapper to call the function without documentation.

518  {
519  // Create a dummy doc_info object
520  DocInfo doc_info;
521  doc_info.directory() = "";
522  doc_info.disable_doc();
523  // Call the other adapt mesh
524  adapt_mesh(doc_info);
525  }
oomph::DocInfo::directory
std::string directory() const
Output directory.
Definition: oomph_utilities.h:524

References oomph::DocInfo::directory(), oomph::DocInfo::disable_doc(), and oomph::RefineableMeshBase::doc_info().

Referenced by adapt(), refine_base_mesh(), refine_selected_elements(), and refine_uniformly().

◆ adapt_mesh() [2/2]

void oomph::TreeBasedRefineableMeshBase::adapt_mesh ( DocInfo &  doc_info)
virtual

Perform the actual tree-based mesh adaptation, documenting the progress in the directory specified in DocInfo object.

Adapt mesh, which exists in two representations, namely as:

  • a FE mesh
  • a forest of Oc or QuadTrees

Refinement/derefinement process is documented (in tecplot-able form) if requested.

Procedure:

  • Loop over all elements and do the refinement for those who want to be refined. Note: Refinement/splitting only allocates elements but doesn't build them.
  • Build the new elements (i.e. give them nodes (create new ones where necessary), assign boundary conditions, and add nodes to mesh and mesh boundaries.
  • For all nodes that were hanging on the previous mesh (and are still marked as such), fill in their nodal values (consistent with the current hanging node scheme) to make sure they are fully functional, should they have become non-hanging during the mesh-adaptation. Then mark the nodes as non-hanging.
  • Unrefine selected elements (which may cause nodes to be re-built).
  • Add the new elements to the mesh (by completely overwriting the old Vector of elements).
  • Delete any nodes that have become obsolete.
  • Mark up hanging nodes and setup hanging node scheme (incl. recursive cleanup for hanging nodes that depend on other hanging nodes).
  • Adjust position of hanging nodes to make sure their position is consistent with the FE-based represenetation of their larger neighbours.
  • run a quick self-test on the neighbour finding scheme and check the integrity of the elements (if PARANOID)
  • doc hanging node status, boundary conditions, neighbour scheme if requested.

After adaptation, all nodes (whether new or old) have up-to-date current and previous values.

If refinement process is being documented, the following information is documented:

  • The files
    • "neighbours.dat"
    • "all_nodes.dat"
    • "new_nodes.dat"
    • "hang_nodes_*.dat" where the * denotes a direction (n,s,e,w) in 2D or (r,l,u,d,f,b) in 3D

can be viewed with

  • QHangingNodes.mcr
  • The file
    • "hangnodes_withmasters.dat"

can be viewed with

  • QHangingNodesWithMasters.mcr

to check the hanging node status.

  • The neighbour status of the elements is documented in
    • "neighbours.dat"

and can be viewed with

  • QuadTreeNeighbours.mcr

Update the boundary element info – this can be a costly procedure and for this reason the mesh writer might have decided not to set up this scheme. If so, we won't change this and suppress its creation...

927  {
928 #ifdef OOMPH_HAS_MPI
929  // Delete any external element storage before performing the adaptation
930  // (in particular, external halo nodes that are on mesh boundaries)
931  this->delete_all_external_storage();
932 #endif
933 
934  // Only perform the adapt step if the mesh has any elements. This is
935  // relevant in a distributed problem with multiple meshes, where a
936  // particular process may not have any elements on a particular submesh.
937  if (this->nelement() > 0)
938  {
939  // Pointer to mesh needs to be passed to some functions
940  Mesh* mesh_pt = this;
941 
942  double t_start = 0.0;
943  if (Global_timings::Doc_comprehensive_timings)
944  {
945  t_start = TimingHelpers::timer();
946  }
947 
948  // Do refinement(=splitting) of elements that have been selected
949  // This function encapsulates the template parameter
950  this->split_elements_if_required();
951 
952 
953  if (Global_timings::Doc_comprehensive_timings)
954  {
955  double t_end = TimingHelpers::timer();
956  oomph_info << "Time for split_elements_if_required: " << t_end - t_start
957  << std::endl;
958  t_start = TimingHelpers::timer();
959  }
960 
961  // Now elements have been created -- build all the leaves
962  //-------------------------------------------------------
963  // Firstly put all the elements into a vector
964  Vector<Tree*> leaf_nodes_pt;
965  Forest_pt->stick_leaves_into_vector(leaf_nodes_pt);
966 
967 
968  if (Global_timings::Doc_comprehensive_timings)
969  {
970  double t_end = TimingHelpers::timer();
971  oomph_info << "Time for stick_leaves_into_vector: " << t_end - t_start
972  << std::endl;
973  t_start = TimingHelpers::timer();
974  }
975 
976  // If we are documenting the output, create the filename
977  std::ostringstream fullname;
978  std::ofstream new_nodes_file;
979  if (doc_info.is_doc_enabled())
980  {
981  fullname << doc_info.directory() << "/new_nodes" << doc_info.number()
982  << ".dat";
983  new_nodes_file.open(fullname.str().c_str());
984  }
985 
986 
987  // Build all elements and store vector of pointers to new nodes
988  // (Note: build() checks if the element has been built
989  // already, i.e. if it's not a new element).
990  Vector<Node*> new_node_pt;
991  bool was_already_built;
992  unsigned long num_tree_nodes = leaf_nodes_pt.size();
993  for (unsigned long e = 0; e < num_tree_nodes; e++)
994  {
995  // Pre-build must be performed before any elements are built
996  leaf_nodes_pt[e]->object_pt()->pre_build(mesh_pt, new_node_pt);
997  }
998  for (unsigned long e = 0; e < num_tree_nodes; e++)
999  {
1000  // Now do the actual build of the new elements
1001  leaf_nodes_pt[e]->object_pt()->build(
1002  mesh_pt, new_node_pt, was_already_built, new_nodes_file);
1003  }
1004 
1005 
1006  double t_end = 0.0;
1007  if (Global_timings::Doc_comprehensive_timings)
1008  {
1009  t_end = TimingHelpers::timer();
1010  oomph_info << "Time for building " << num_tree_nodes
1011  << " new elements: " << t_end - t_start << std::endl;
1012  t_start = TimingHelpers::timer();
1013  }
1014 
1015  // Close the new nodes files, if it was opened
1016  if (doc_info.is_doc_enabled())
1017  {
1018  new_nodes_file.close();
1019  }
1020 
1021  // Loop over all nodes in mesh and free the dofs of those that were
1022  //-----------------------------------------------------------------
1023  // pinned only because they were hanging nodes. Also update their
1024  //-----------------------------------------------------------------
1025  // nodal values so that they contain data that is consistent
1026  //----------------------------------------------------------
1027  // with the hanging node representation
1028  //-------------------------------------
1029  // (Even if the nodal data isn't actually accessed because the node
1030  // is still hanging -- we don't know this yet, and this step makes
1031  // sure that all nodes are fully functional and up-to-date, should
1032  // they become non-hanging below).
1033  //
1034  //
1035  // However, if we have a fixed mesh and hanging nodes on the boundary
1036  // become non-hanging they will not necessarily respect the curvilinear
1037  // boundaries. This can only happen in 3D of course because it is not
1038  // possible to have hanging nodes on boundaries in 2D.
1039  // The solution is to store those nodes on the boundaries that are
1040  // currently hanging and then check to see whether they have changed
1041  // status at the end of the refinement procedure.
1042  // If it has changed, then we need to adjust their positions.
1043  const unsigned n_boundary = this->nboundary();
1044  const unsigned mesh_dim = this->finite_element_pt(0)->dim();
1045  Vector<std::set<Node*>> hanging_nodes_on_boundary_pt(n_boundary);
1046 
1047  unsigned long n_node = this->nnode();
1048  for (unsigned long n = 0; n < n_node; n++)
1049  {
1050  // Get the pointer to the node
1051  Node* nod_pt = this->node_pt(n);
1052 
1053  // Get the number of values in the node
1054  unsigned n_value = nod_pt->nvalue();
1055 
1056  // We need to find if any of the values are hanging
1057  bool is_hanging = nod_pt->is_hanging();
1058  // Loop over the values and find out whether any are hanging
1059  for (unsigned n = 0; n < n_value; n++)
1060  {
1061  is_hanging |= nod_pt->is_hanging(n);
1062  }
1063 
1064  // If the node is hanging then ...
1065  if (is_hanging)
1066  {
1067  // Unless they are turned into hanging nodes again below
1068  // (this might or might not happen), fill in all the necessary
1069  // data to make them 'proper' nodes again.
1070 
1071  // Reconstruct the nodal values/position from the node's
1072  // hanging node representation
1073  unsigned nt = nod_pt->ntstorage();
1074  Vector<double> values(n_value);
1075  unsigned n_dim = nod_pt->ndim();
1076  Vector<double> position(n_dim);
1077  // Loop over all history values
1078  for (unsigned t = 0; t < nt; t++)
1079  {
1080  nod_pt->value(t, values);
1081  for (unsigned i = 0; i < n_value; i++)
1082  {
1083  nod_pt->set_value(t, i, values[i]);
1084  }
1085  nod_pt->position(t, position);
1086  for (unsigned i = 0; i < n_dim; i++)
1087  {
1088  nod_pt->x(t, i) = position[i];
1089  }
1090  }
1091 
1092  // If it's an algebraic node: Update its previous nodal positions too
1093  AlgebraicNode* alg_node_pt = dynamic_cast<AlgebraicNode*>(nod_pt);
1094  if (alg_node_pt != 0)
1095  {
1096  bool update_all_time_levels = true;
1097  alg_node_pt->node_update(update_all_time_levels);
1098  }
1099 
1100 
1101  // If it's a Solid node, update Lagrangian coordinates
1102  // from its hanging node representation
1103  SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
1104  if (solid_node_pt != 0)
1105  {
1106  unsigned n_lagrangian = solid_node_pt->nlagrangian();
1107  for (unsigned i = 0; i < n_lagrangian; i++)
1108  {
1109  solid_node_pt->xi(i) = solid_node_pt->lagrangian_position(i);
1110  }
1111  }
1112 
1113  // Now store geometrically hanging nodes on boundaries that
1114  // may need updating after refinement.
1115  // There will only be a problem if we have 3 spatial dimensions
1116  if ((mesh_dim > 2) && (nod_pt->is_hanging()))
1117  {
1118  // If the node is on a boundary then add a pointer to the node
1119  // to our lookup scheme
1120  if (nod_pt->is_on_boundary())
1121  {
1122  // Storage for the boundaries on which the Node is located
1123  std::set<unsigned>* boundaries_pt;
1124  nod_pt->get_boundaries_pt(boundaries_pt);
1125  if (boundaries_pt != 0)
1126  {
1127  // Loop over the boundaries and add a pointer to the node
1128  // to the appropriate storage scheme
1129  for (std::set<unsigned>::iterator it = boundaries_pt->begin();
1130  it != boundaries_pt->end();
1131  ++it)
1132  {
1133  hanging_nodes_on_boundary_pt[*it].insert(nod_pt);
1134  }
1135  }
1136  }
1137  }
1138 
1139  } // End of is_hanging
1140 
1141  // Initially mark all nodes as 'non-hanging' and `obsolete'
1142  nod_pt->set_nonhanging();
1143  nod_pt->set_obsolete();
1144  }
1145 
1146  if (Global_timings::Doc_comprehensive_timings)
1147  {
1148  t_end = TimingHelpers::timer();
1149  oomph_info << "Time for sorting out initial hanging status: "
1150  << t_end - t_start << std::endl;
1151  t_start = TimingHelpers::timer();
1152  }
1153 
1154  // Unrefine all the selected elements: This needs to be
1155  //-----------------------------------------------------
1156  // all elements, because the father elements are not actually leaves.
1157  //-------------------------------------------------------------------
1158 
1159  // Unrefine
1160  for (unsigned long e = 0; e < Forest_pt->ntree(); e++)
1161  {
1162  Forest_pt->tree_pt(e)->traverse_all(&Tree::merge_sons_if_required,
1163  mesh_pt);
1164  }
1165 
1166  if (Global_timings::Doc_comprehensive_timings)
1167  {
1168  t_end = TimingHelpers::timer();
1169  oomph_info << "Time for unrefinement: " << t_end - t_start << std::endl;
1170  t_start = TimingHelpers::timer();
1171  }
1172 
1173  // Add the newly created elements to mesh
1174  //---------------------------------------
1175 
1176  // Stick all elements into a new vector
1177  //(note the leaves may have changed, so this is not duplicated work)
1178  Vector<Tree*> tree_nodes_pt;
1179  Forest_pt->stick_leaves_into_vector(tree_nodes_pt);
1180 
1181  // Copy the elements into the mesh Vector
1182  num_tree_nodes = tree_nodes_pt.size();
1183  Element_pt.resize(num_tree_nodes);
1184  for (unsigned long e = 0; e < num_tree_nodes; e++)
1185  {
1186  Element_pt[e] = tree_nodes_pt[e]->object_pt();
1187 
1188  // Now loop over all nodes in element and mark them as non-obsolete
1189  // Logic: Initially all nodes in the unrefined mesh were labeled
1190  // as deleteable. Then we create new elements (whose newly created
1191  // nodes are obviously non-obsolete), and killed some other elements (by
1192  // by deleting them and marking the nodes that were not shared by
1193  // their father as obsolete. Now we loop over all the remaining
1194  // elements and (re-)label all their nodes as non-obsolete. This
1195  // saves some nodes that were regarded as obsolete by deleted
1196  // elements but are still required in some surviving ones
1197  // from a tragic early death...
1198  FiniteElement* this_el_pt = this->finite_element_pt(e);
1199  unsigned n_node = this_el_pt->nnode(); // caching pre-loop
1200  for (unsigned n = 0; n < n_node; n++)
1201  {
1202  this_el_pt->node_pt(n)->set_non_obsolete();
1203  }
1204  }
1205 
1206 
1207  if (Global_timings::Doc_comprehensive_timings)
1208  {
1209  t_end = TimingHelpers::timer();
1210  oomph_info << "Time for adding elements to mesh: " << t_end - t_start
1211  << std::endl;
1212  t_start = TimingHelpers::timer();
1213  }
1214 
1215  // Cannot delete nodes that are still marked as obsolete
1216  // because they may still be required to assemble the hanging schemes
1217  //-------------------------------------------------------------------
1218 
1219  // Mark up hanging nodes
1220  //----------------------
1221 
1222  // Output streams for the hanging nodes
1223  Vector<std::ofstream*> hanging_output_files;
1224  // Setup the output files for hanging nodes, this must be called
1225  // precisely once for the forest. Note that the files will only
1226  // actually be opened if doc_info.is_doc_enabled() is true
1227  Forest_pt->open_hanging_node_files(doc_info, hanging_output_files);
1228 
1229  for (unsigned long e = 0; e < num_tree_nodes; e++)
1230  {
1231  // Generic setup
1232  tree_nodes_pt[e]->object_pt()->setup_hanging_nodes(
1233  hanging_output_files);
1234  // Element specific setup
1235  tree_nodes_pt[e]->object_pt()->further_setup_hanging_nodes();
1236  }
1237 
1238  // Close the hanging node files and delete the memory allocated
1239  // for the streams
1240  Forest_pt->close_hanging_node_files(doc_info, hanging_output_files);
1241 
1242 
1243  if (Global_timings::Doc_comprehensive_timings)
1244  {
1245  t_end = TimingHelpers::timer();
1246  oomph_info << "Time for setup_hanging_nodes() and "
1247  "further_setup_hanging_nodes() for "
1248  << num_tree_nodes << " elements: " << t_end - t_start
1249  << std::endl;
1250  t_start = TimingHelpers::timer();
1251  }
1252 
1253  // Read out the number of continously interpolated values
1254  // from one of the elements (assuming it's the same in all elements)
1255  unsigned ncont_interpolated_values =
1256  tree_nodes_pt[0]->object_pt()->ncont_interpolated_values();
1257 
1258  // Complete the hanging nodes schemes by dealing with the
1259  // recursively hanging nodes
1260  complete_hanging_nodes(ncont_interpolated_values);
1261 
1262  if (Global_timings::Doc_comprehensive_timings)
1263  {
1264  t_end = TimingHelpers::timer();
1265  oomph_info << "Time for complete_hanging_nodes: " << t_end - t_start
1266  << std::endl;
1267  t_start = TimingHelpers::timer();
1268  }
1269 
1274  if (Lookup_for_elements_next_boundary_is_setup)
1275  {
1276  this->setup_boundary_element_info();
1277  }
1278 
1279  if (Global_timings::Doc_comprehensive_timings)
1280  {
1281  t_end = TimingHelpers::timer();
1282  oomph_info << "Time for boundary element info: " << t_end - t_start
1283  << std::endl;
1284  t_start = TimingHelpers::timer();
1285  }
1286 
1287 #ifdef PARANOID
1288 
1289  // Doc/check the neighbours
1290  //-------------------------
1291  Vector<Tree*> all_tree_nodes_pt;
1292  Forest_pt->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
1293 
1294  // Check the neighbours
1295  Forest_pt->check_all_neighbours(doc_info);
1296 
1297  // Check the integrity of the elements
1298  // -----------------------------------
1299 
1300  // Loop over elements and get the elemental integrity
1301  double max_error = 0.0;
1302  for (unsigned long e = 0; e < num_tree_nodes; e++)
1303  {
1304  double max_el_error;
1305  tree_nodes_pt[e]->object_pt()->check_integrity(max_el_error);
1306  // If the elemental error is greater than our maximum error
1307  // reset the maximum
1308  if (max_el_error > max_error)
1309  {
1310  max_error = max_el_error;
1311  }
1312  }
1313 
1314  if (max_error > RefineableElement::max_integrity_tolerance())
1315  {
1316  std::ostringstream error_stream;
1317  error_stream << "Mesh refined: Max. error in integrity check: "
1318  << max_error << " is too big\n";
1319  error_stream
1320  << "i.e. bigger than RefineableElement::max_integrity_tolerance()="
1321  << RefineableElement::max_integrity_tolerance() << "\n"
1322  << std::endl;
1323 
1324  std::ofstream some_file;
1325  some_file.open("ProblemMesh.dat");
1326  for (unsigned long n = 0; n < n_node; n++)
1327  {
1328  // Get the pointer to the node
1329  Node* nod_pt = this->node_pt(n);
1330  // Get the dimension
1331  unsigned n_dim = nod_pt->ndim();
1332  // Output the coordinates
1333  for (unsigned i = 0; i < n_dim; i++)
1334  {
1335  some_file << this->node_pt(n)->x(i) << " ";
1336  }
1337  some_file << std::endl;
1338  }
1339  some_file.close();
1340 
1341  error_stream << "Doced problem mesh in ProblemMesh.dat" << std::endl;
1342 
1343  throw OomphLibError(
1344  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1345  }
1346  else
1347  {
1348  oomph_info << "Mesh refined: Max. error in integrity check: "
1349  << max_error << " is OK" << std::endl;
1350  oomph_info
1351  << "i.e. less than RefineableElement::max_integrity_tolerance()="
1352  << RefineableElement::max_integrity_tolerance() << "\n"
1353  << std::endl;
1354  }
1355 
1356 
1357  if (Global_timings::Doc_comprehensive_timings)
1358  {
1359  t_end = TimingHelpers::timer();
1360  oomph_info << "Time for (paranoid only) checking of integrity: "
1361  << t_end - t_start << std::endl;
1362  t_start = TimingHelpers::timer();
1363  }
1364 
1365 #endif
1366 
1367  // Loop over all elements other than the final level and deactivate the
1368  // objects, essentially set the pointer that point to nodes that are
1369  // about to be deleted to NULL. This must take place here because nodes
1370  // addressed by elements that are dead but still living in the tree might
1371  // have been made obsolete in the last round of refinement
1372  for (unsigned long e = 0; e < Forest_pt->ntree(); e++)
1373  {
1374  Forest_pt->tree_pt(e)->traverse_all_but_leaves(
1375  &Tree::deactivate_object);
1376  }
1377 
1378  // Now we can prune the dead nodes from the mesh.
1379  Vector<Node*> deleted_node_pt = this->prune_dead_nodes();
1380 
1381  if (Global_timings::Doc_comprehensive_timings)
1382  {
1383  t_end = TimingHelpers::timer();
1384  oomph_info << "Time for deactivating objects and pruning nodes: "
1385  << t_end - t_start << std::endl;
1386  t_start = TimingHelpers::timer();
1387  }
1388 
1389  // Finally: Reorder the nodes within the mesh's node vector
1390  // to establish a standard ordering regardless of the sequence
1391  // of mesh refinements -- this is required to allow dump/restart
1392  // on refined meshes
1393  this->reorder_nodes();
1394 
1395  if (Global_timings::Doc_comprehensive_timings)
1396  {
1397  t_end = TimingHelpers::timer();
1398  oomph_info << "Time for reordering " << nnode()
1399  << " nodes: " << t_end - t_start << std::endl;
1400  t_start = TimingHelpers::timer();
1401  }
1402 
1403  // Now we can correct the nodes on boundaries that were hanging that
1404  // are no longer hanging
1405  // Only bother if we have more than two dimensions
1406  if (mesh_dim > 2)
1407  {
1408  // Loop over the boundaries
1409  for (unsigned b = 0; b < n_boundary; b++)
1410  {
1411  // Remove deleted nodes from the set
1412  unsigned n_del = deleted_node_pt.size();
1413  for (unsigned j = 0; j < n_del; j++)
1414  {
1415  hanging_nodes_on_boundary_pt[b].erase(deleted_node_pt[j]);
1416  }
1417 
1418  // If the nodes that were hanging are still hanging then remove them
1419  // from the set (note increment is not in for command for efficiencty)
1420  for (std::set<Node*>::iterator it =
1421  hanging_nodes_on_boundary_pt[b].begin();
1422  it != hanging_nodes_on_boundary_pt[b].end();)
1423  {
1424  if ((*it)->is_hanging())
1425  {
1426  hanging_nodes_on_boundary_pt[b].erase(it++);
1427  }
1428  else
1429  {
1430  ++it;
1431  }
1432  }
1433 
1434  // Are there any nodes that have changed geometric hanging status
1435  // on the boundary
1436  // The slightly painful part is that we must adjust the position
1437  // via the macro-elements which are only available through the
1438  // elements and not the nodes.
1439  if (hanging_nodes_on_boundary_pt[b].size() > 0)
1440  {
1441  // If so we loop over all elements adjacent to the boundary
1442  unsigned n_boundary_element = this->nboundary_element(b);
1443  for (unsigned e = 0; e < n_boundary_element; ++e)
1444  {
1445  // Get a pointer to the element
1446  FiniteElement* el_pt = this->boundary_element_pt(b, e);
1447 
1448  // Do we have a solid element
1449  SolidFiniteElement* solid_el_pt =
1450  dynamic_cast<SolidFiniteElement*>(el_pt);
1451 
1452  // Determine whether there is a macro element
1453  bool macro_present = (el_pt->macro_elem_pt() != 0);
1454  // Or a solid macro element
1455  if (solid_el_pt != 0)
1456  {
1457  macro_present |= (solid_el_pt->undeformed_macro_elem_pt() != 0);
1458  }
1459 
1460  // Only bother to do anything if there is a macro element
1461  // or undeformed macro element in a SolidElement
1462  if (macro_present)
1463  {
1464  // Loop over the nodes
1465  // ALH: (could optimise to only loop over
1466  // node associated with the boundary with more effort)
1467  unsigned n_el_node = el_pt->nnode();
1468  for (unsigned n = 0; n < n_el_node; n++)
1469  {
1470  // Cache pointer to the node
1471  Node* nod_pt = el_pt->node_pt(n);
1472  if (nod_pt->is_on_boundary(b))
1473  {
1474  // Is the Node in our set
1475  std::set<Node*>::iterator it =
1476  hanging_nodes_on_boundary_pt[b].find(nod_pt);
1477 
1478  // If we have found the Node then update the position
1479  // to be consistent with the macro-element representation
1480  if (it != hanging_nodes_on_boundary_pt[b].end())
1481  {
1482  // Specialise local and global coordinates to 3D
1483  // because there is only a problem in 3D.
1484  Vector<double> s(3), x(3);
1485 
1486  // Find the local coordinate of the ndoe
1487  el_pt->local_coordinate_of_node(n, s);
1488 
1489  // Find the number of time history values
1490  const unsigned ntstorage = nod_pt->ntstorage();
1491 
1492  // Do we have a solid node
1493  SolidNode* solid_node_pt =
1494  dynamic_cast<SolidNode*>(nod_pt);
1495  if (solid_node_pt)
1496  {
1497  // Assign Lagrangian coordinates from undeformed
1498  // macro element (if it has one -- get_x_and_xi()
1499  // does "the right thing" anyway. Leave actual
1500  // nodal positions alone -- we're doing a solid
1501  // mechanics problem and once we're going
1502  // the nodal positions are always computed, never
1503  // (automatically) reset to macro-element based
1504  // positions; not even on pinned boundaries
1505  // because the user may have other ideas about where
1506  // these should go -- pinning means "don't touch the
1507  // value", not "leave where the macro-element thinks
1508  // it should be"
1509  Vector<double> x_fe(3), xi(3), xi_fe(3);
1510  solid_el_pt->get_x_and_xi(s, x_fe, x, xi_fe, xi);
1511  for (unsigned i = 0; i < 3; i++)
1512  {
1513  solid_node_pt->xi(i) = xi[i];
1514  }
1515  }
1516  else
1517  {
1518  // Set position and history values from the
1519  // macro-element representation
1520  for (unsigned t = 0; t < ntstorage; t++)
1521  {
1522  // Get the history value from the macro element
1523  el_pt->get_x(t, s, x);
1524 
1525  // Set the coordinate to that of the macroelement
1526  // representation
1527  for (unsigned i = 0; i < 3; i++)
1528  {
1529  nod_pt->x(t, i) = x[i];
1530  }
1531  }
1532  } // End of non-solid node case
1533 
1534  // Now remove the node from the list
1535  hanging_nodes_on_boundary_pt[b].erase(it);
1536 
1537  // If there are no Nodes left then exit the loops
1538  if (hanging_nodes_on_boundary_pt[b].size() == 0)
1539  {
1540  e = n_boundary_element;
1541  break;
1542  }
1543  }
1544  }
1545  }
1546  } // End of macro element case
1547  }
1548  }
1549  }
1550  } // End of case when we have fixed nodal positions
1551 
1552 
1553  // Final doc
1554  //-----------
1555  if (doc_info.is_doc_enabled())
1556  {
1557  // Doc the boundary conditions ('0' for non-existent, '1' for free,
1558  //----------------------------------------------------------------
1559  // '2' for pinned -- ideal for tecplot scatter sizing.
1560  //----------------------------------------------------
1561  // num_tree_nodes=tree_nodes_pt.size();
1562 
1563  // Determine maximum number of values at any node in this type of
1564  // element
1565  RefineableElement* el_pt = tree_nodes_pt[0]->object_pt();
1566  // Initalise max_nval
1567  unsigned max_nval = 0;
1568  for (unsigned n = 0; n < el_pt->nnode(); n++)
1569  {
1570  if (el_pt->node_pt(n)->nvalue() > max_nval)
1571  {
1572  max_nval = el_pt->node_pt(n)->nvalue();
1573  }
1574  }
1575 
1576  // Open the output file
1577  std::ofstream bcs_file;
1578  fullname.str("");
1579  fullname << doc_info.directory() << "/bcs" << doc_info.number()
1580  << ".dat";
1581  bcs_file.open(fullname.str().c_str());
1582 
1583  // Loop over elements
1584  for (unsigned long e = 0; e < num_tree_nodes; e++)
1585  {
1586  el_pt = tree_nodes_pt[e]->object_pt();
1587  // Loop over nodes in element
1588  unsigned n_nod = el_pt->nnode();
1589  for (unsigned n = 0; n < n_nod; n++)
1590  {
1591  // Get pointer to the node
1592  Node* nod_pt = el_pt->node_pt(n);
1593  // Find the dimension of the node
1594  unsigned n_dim = nod_pt->ndim();
1595  // Write the nodal coordinates to the file
1596  for (unsigned i = 0; i < n_dim; i++)
1597  {
1598  bcs_file << nod_pt->x(i) << " ";
1599  }
1600 
1601  // Loop over all values in this element
1602  for (unsigned i = 0; i < max_nval; i++)
1603  {
1604  // Value exists at this node:
1605  if (i < nod_pt->nvalue())
1606  {
1607  bcs_file << " " << 1 + nod_pt->is_pinned(i);
1608  }
1609  // ...if not just dump out a zero
1610  else
1611  {
1612  bcs_file << " 0 ";
1613  }
1614  }
1615  bcs_file << std::endl;
1616  }
1617  }
1618  bcs_file.close();
1619 
1620  // Doc all nodes
1621  //---------------
1622  std::ofstream all_nodes_file;
1623  fullname.str("");
1624  fullname << doc_info.directory() << "/all_nodes" << doc_info.number()
1625  << ".dat";
1626  all_nodes_file.open(fullname.str().c_str());
1627 
1628  all_nodes_file << "ZONE \n";
1629 
1630  // Need to recompute the number of nodes since it may have
1631  // changed during mesh refinement/unrefinement
1632  n_node = this->nnode();
1633  for (unsigned long n = 0; n < n_node; n++)
1634  {
1635  Node* nod_pt = this->node_pt(n);
1636  unsigned n_dim = nod_pt->ndim();
1637  for (unsigned i = 0; i < n_dim; i++)
1638  {
1639  all_nodes_file << this->node_pt(n)->x(i) << " ";
1640  }
1641  all_nodes_file << std::endl;
1642  }
1643 
1644  all_nodes_file.close();
1645 
1646 
1647  // Doc all hanging nodes:
1648  //-----------------------
1649  std::ofstream some_file;
1650  fullname.str("");
1651  fullname << doc_info.directory() << "/all_hangnodes"
1652  << doc_info.number() << ".dat";
1653  some_file.open(fullname.str().c_str());
1654  for (unsigned long n = 0; n < n_node; n++)
1655  {
1656  Node* nod_pt = this->node_pt(n);
1657 
1658  if (nod_pt->is_hanging())
1659  {
1660  unsigned n_dim = nod_pt->ndim();
1661  for (unsigned i = 0; i < n_dim; i++)
1662  {
1663  some_file << nod_pt->x(i) << " ";
1664  }
1665 
1666  // ALH: Added this to stop Solid problems seg-faulting
1667  if (this->node_pt(n)->nvalue() > 0)
1668  {
1669  some_file << " " << nod_pt->raw_value(0);
1670  }
1671  some_file << std::endl;
1672  }
1673  }
1674  some_file.close();
1675 
1676  // Doc all hanging nodes and their masters
1677  // View with QHangingNodesWithMasters.mcr
1678  fullname.str("");
1679  fullname << doc_info.directory() << "/geometric_hangnodes_withmasters"
1680  << doc_info.number() << ".dat";
1681  some_file.open(fullname.str().c_str());
1682  for (unsigned long n = 0; n < n_node; n++)
1683  {
1684  Node* nod_pt = this->node_pt(n);
1685  if (nod_pt->is_hanging())
1686  {
1687  unsigned n_dim = nod_pt->ndim();
1688  unsigned nmaster = nod_pt->hanging_pt()->nmaster();
1689  some_file << "ZONE I=" << nmaster + 1 << std::endl;
1690  for (unsigned i = 0; i < n_dim; i++)
1691  {
1692  some_file << nod_pt->x(i) << " ";
1693  }
1694  some_file << " 2 " << std::endl;
1695 
1696  for (unsigned imaster = 0; imaster < nmaster; imaster++)
1697  {
1698  Node* master_nod_pt =
1699  nod_pt->hanging_pt()->master_node_pt(imaster);
1700  unsigned n_dim = master_nod_pt->ndim();
1701  for (unsigned i = 0; i < n_dim; i++)
1702  {
1703  some_file << master_nod_pt->x(i) << " ";
1704  }
1705  some_file << " 1 " << std::endl;
1706  }
1707  }
1708  }
1709  some_file.close();
1710 
1711  // Doc all hanging nodes and their masters
1712  // View with QHangingNodesWithMasters.mcr
1713  for (unsigned i = 0; i < ncont_interpolated_values; i++)
1714  {
1715  fullname.str("");
1716  fullname << doc_info.directory()
1717  << "/nonstandard_hangnodes_withmasters" << i << "_"
1718  << doc_info.number() << ".dat";
1719  some_file.open(fullname.str().c_str());
1720  unsigned n_nod = this->nnode();
1721  for (unsigned long n = 0; n < n_nod; n++)
1722  {
1723  Node* nod_pt = this->node_pt(n);
1724  if (nod_pt->is_hanging(i))
1725  {
1726  if (nod_pt->hanging_pt(i) != nod_pt->hanging_pt())
1727  {
1728  unsigned nmaster = nod_pt->hanging_pt(i)->nmaster();
1729  some_file << "ZONE I=" << nmaster + 1 << std::endl;
1730  unsigned n_dim = nod_pt->ndim();
1731  for (unsigned j = 0; j < n_dim; j++)
1732  {
1733  some_file << nod_pt->x(j) << " ";
1734  }
1735  some_file << " 2 " << std::endl;
1736  for (unsigned imaster = 0; imaster < nmaster; imaster++)
1737  {
1738  Node* master_nod_pt =
1739  nod_pt->hanging_pt(i)->master_node_pt(imaster);
1740  unsigned n_dim = master_nod_pt->ndim();
1741  for (unsigned j = 0; j < n_dim; j++)
1742  {
1743  // some_file << master_nod_pt->x(i) << " ";
1744  }
1745  some_file << " 1 " << std::endl;
1746  }
1747  }
1748  }
1749  }
1750  some_file.close();
1751  }
1752  } // End of documentation
1753  } // End if (this->nelement()>0)
1754 
1755 
1756 #ifdef OOMPH_HAS_MPI
1757 
1758  // Now (re-)classify halo and haloed nodes and synchronise hanging
1759  // nodes
1760  if (this->is_mesh_distributed())
1761  {
1762  classify_halo_and_haloed_nodes(doc_info, doc_info.is_doc_enabled());
1763  }
1764 
1765 #endif
1766  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
b
Scalar * b
Definition: benchVecAdd.cpp:17
oomph::DocInfo::number
unsigned & number()
Number used (e.g.) for labeling output files.
Definition: oomph_utilities.h:554
oomph::FiniteElement::dim
unsigned dim() const
Definition: elements.h:2611
oomph::Mesh::Lookup_for_elements_next_boundary_is_setup
bool Lookup_for_elements_next_boundary_is_setup
Definition: mesh.h:87
oomph::Mesh::finite_element_pt
FiniteElement * finite_element_pt(const unsigned &e) const
Definition: mesh.h:473
oomph::Mesh::nboundary_element
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
oomph::Mesh::setup_boundary_element_info
virtual void setup_boundary_element_info()
Definition: mesh.h:275
oomph::Mesh::nboundary
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
oomph::Mesh::boundary_element_pt
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840
oomph::Mesh::nnode
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
oomph::Mesh::Mesh
Mesh()
Default constructor.
Definition: mesh.h:236
oomph::Mesh::node_pt
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
oomph::Mesh::prune_dead_nodes
Vector< Node * > prune_dead_nodes()
Definition: mesh.cc:966
oomph::Mesh::Element_pt
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
oomph::Node::x
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
oomph::RefineableElement::max_integrity_tolerance
static double & max_integrity_tolerance()
Max. allowed discrepancy in element integrity check.
Definition: refineable_elements.h:495
oomph::RefineableMeshBase::max_error
double & max_error()
Definition: refineable_mesh.h:177
oomph::TreeBasedRefineableMeshBase::split_elements_if_required
virtual void split_elements_if_required()=0
oomph::TreeBasedRefineableMeshBase::complete_hanging_nodes
void complete_hanging_nodes(const int &ncont_interpolated_values)
Complete the hanging node scheme recursively.
Definition: refineable_mesh.cc:2271
oomph::TreeForest::check_all_neighbours
virtual void check_all_neighbours(DocInfo &doc_info)=0
oomph::TreeForest::stick_all_tree_nodes_into_vector
void stick_all_tree_nodes_into_vector(Vector< Tree * > &all_forest_nodes)
Traverse forest and stick pointers to all "nodes" into Vector.
Definition: tree.cc:405
oomph::TreeForest::open_hanging_node_files
virtual void open_hanging_node_files(DocInfo &doc_info, Vector< std::ofstream * > &output_stream)=0
oomph::TreeForest::ntree
unsigned ntree()
Number of trees in forest.
Definition: tree.h:460
oomph::TreeForest::stick_leaves_into_vector
void stick_leaves_into_vector(Vector< Tree * > &forest_nodes)
Traverse forst and stick pointers to leaf "nodes" into Vector.
Definition: tree.cc:379
oomph::TreeForest::tree_pt
TreeRoot * tree_pt(const unsigned &i) const
Return pointer to i-th tree in forest.
Definition: tree.h:466
oomph::TreeForest::close_hanging_node_files
void close_hanging_node_files(DocInfo &doc_info, Vector< std::ofstream * > &output_stream)
Definition: tree.cc:432
oomph::Tree::traverse_all
void traverse_all(Tree::VoidMemberFctPt member_function)
Definition: tree.cc:145
oomph::Tree::deactivate_object
void deactivate_object()
Call the RefineableElement's deactivate_element() function.
Definition: tree.cc:341
oomph::Tree::traverse_all_but_leaves
void traverse_all_but_leaves(Tree::VoidMemberFctPt member_function)
Definition: tree.cc:184
oomph::Tree::merge_sons_if_required
void merge_sons_if_required(Mesh *&mesh_pt)
Definition: tree.cc:302
Eigen::placeholders::end
static constexpr lastp1_t end
Definition: IndexedViewHelper.h:79
s
RealScalar s
Definition: level1_cplx_impl.h:130
oomph::Global_timings::Doc_comprehensive_timings
bool Doc_comprehensive_timings
Definition: oomph_definitions.cc:49
oomph::TimingHelpers::timer
double timer()
returns the time in seconds after some point in past
Definition: oomph_utilities.cc:1295
plotDoE.x
list x
Definition: plotDoE.py:28
plotPSD.t
t
Definition: plotPSD.py:36
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References b, oomph::Mesh::boundary_element_pt(), oomph::TreeForest::check_all_neighbours(), oomph::TreeForest::close_hanging_node_files(), complete_hanging_nodes(), oomph::Tree::deactivate_object(), oomph::Mesh::delete_all_external_storage(), oomph::FiniteElement::dim(), oomph::DocInfo::directory(), oomph::Global_timings::Doc_comprehensive_timings, oomph::RefineableMeshBase::doc_info(), e(), oomph::Mesh::Element_pt, Eigen::placeholders::end, oomph::Mesh::finite_element_pt(), Forest_pt, oomph::Node::get_boundaries_pt(), oomph::FiniteElement::get_x(), oomph::SolidFiniteElement::get_x_and_xi(), oomph::Node::hanging_pt(), i, oomph::DocInfo::is_doc_enabled(), oomph::Node::is_hanging(), oomph::Mesh::is_mesh_distributed(), oomph::Node::is_on_boundary(), oomph::Data::is_pinned(), j, oomph::SolidNode::lagrangian_position(), oomph::FiniteElement::local_coordinate_of_node(), oomph::Mesh::Lookup_for_elements_next_boundary_is_setup, oomph::FiniteElement::macro_elem_pt(), oomph::HangInfo::master_node_pt(), oomph::RefineableMeshBase::max_error(), oomph::RefineableElement::max_integrity_tolerance(), oomph::Tree::merge_sons_if_required(), n, oomph::Mesh::nboundary(), oomph::Mesh::nboundary_element(), oomph::Node::ndim(), oomph::Mesh::nelement(), oomph::SolidNode::nlagrangian(), oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::Mesh::nnode(), oomph::FiniteElement::node_pt(), oomph::Mesh::node_pt(), oomph::AlgebraicNode::node_update(), oomph::TreeForest::ntree(), oomph::Data::ntstorage(), oomph::DocInfo::number(), oomph::Data::nvalue(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::TreeForest::open_hanging_node_files(), oomph::Node::position(), oomph::Mesh::prune_dead_nodes(), oomph::Node::raw_value(), oomph::Mesh::reorder_nodes(), s, oomph::Node::set_non_obsolete(), oomph::Node::set_nonhanging(), oomph::Node::set_obsolete(), oomph::Data::set_value(), oomph::Mesh::setup_boundary_element_info(), size, split_elements_if_required(), oomph::TreeForest::stick_all_tree_nodes_into_vector(), oomph::TreeForest::stick_leaves_into_vector(), plotPSD::t, oomph::TimingHelpers::timer(), oomph::Tree::traverse_all(), oomph::Tree::traverse_all_but_leaves(), oomph::TreeForest::tree_pt(), oomph::SolidFiniteElement::undeformed_macro_elem_pt(), oomph::Node::value(), plotDoE::x, oomph::Node::x(), and oomph::SolidNode::xi().

◆ complete_hanging_nodes()

void oomph::TreeBasedRefineableMeshBase::complete_hanging_nodes ( const int &  ncont_interpolated_values)
protected

Complete the hanging node scheme recursively.

Complete the hanging node scheme recursively. After the initial markup scheme, hanging nodes can depend on other hanging nodes —> AAAAAAAAARGH! Need to translate this into a scheme where all hanging nodes only depend on non-hanging nodes...

2273  {
2274  // Number of nodes in mesh
2275  unsigned long n_node = this->nnode();
2276  double min_weight = 1.0e-8; // RefineableBrickElement::min_weight_value();
2277 
2278  // Loop over the nodes in the mesh
2279  for (unsigned long n = 0; n < n_node; n++)
2280  {
2281  // Assign a local pointer to the node
2282  Node* nod_pt = this->node_pt(n);
2283 
2284  // Loop over the values,
2285  // N.B. geometric hanging data is stored at the index -1
2286  for (int i = -1; i < ncont_interpolated_values; i++)
2287  {
2288  // Is the node hanging?
2289  if (nod_pt->is_hanging(i))
2290  {
2291  // If it is geometric OR has hanging node data that differs
2292  // from the geometric data, we must do some work
2293  if ((i == -1) || (nod_pt->hanging_pt(i) != nod_pt->hanging_pt()))
2294  {
2295  // Find out the ultimate map of dependencies: Master nodes
2296  // and associated weights
2297  Vector<Node*> master_nodes;
2298  Vector<double> hanging_weights;
2299  complete_hanging_nodes_recursively(
2300  nod_pt, master_nodes, hanging_weights, i);
2301 
2302  // put them into a map to merge all the occurences of the same node
2303  // (add the weights)
2304  std::map<Node*, double> hang_weights;
2305  unsigned n_master = master_nodes.size();
2306  for (unsigned k = 0; k < n_master; k++)
2307  {
2308  if (std::fabs(hanging_weights[k]) > min_weight)
2309  hang_weights[master_nodes[k]] += hanging_weights[k];
2310  }
2311 
2312  // Create new hanging data (we know how many data there are)
2313  HangInfo* hang_pt = new HangInfo(hang_weights.size());
2314 
2315  unsigned hang_weights_index = 0;
2316  // Copy the map into the HangInfo object
2317  typedef std::map<Node*, double>::iterator IT;
2318  for (IT it = hang_weights.begin(); it != hang_weights.end(); ++it)
2319  {
2320  hang_pt->set_master_node_pt(
2321  hang_weights_index, it->first, it->second);
2322  ++hang_weights_index;
2323  }
2324 
2325  // Assign the new hanging pointer to the appropriate value
2326  nod_pt->set_hanging_pt(hang_pt, i);
2327  }
2328  }
2329  }
2330  }
2331 
2332 #ifdef PARANOID
2333 
2334  // Check hanging node scheme: The weights need to add up to one
2335  //-------------------------------------------------------------
2336  // Loop over all values indices
2337  for (int i = -1; i < ncont_interpolated_values; i++)
2338  {
2339  // Loop over all nodes in mesh
2340  for (unsigned long n = 0; n < n_node; n++)
2341  {
2342  // Set a local pointer to the node
2343  Node* nod_pt = this->node_pt(n);
2344 
2345  // Is it hanging?
2346  if (nod_pt->is_hanging(i))
2347  {
2348  unsigned nmaster = nod_pt->hanging_pt(i)->nmaster();
2349  double sum = 0.0;
2350  for (unsigned imaster = 0; imaster < nmaster; imaster++)
2351  {
2352  sum += nod_pt->hanging_pt(i)->master_weight(imaster);
2353  }
2354  if (std::fabs(sum - 1.0) > 1.0e-7)
2355  {
2356  oomph_info << "WARNING: Sum of master node weights fabs(sum-1.0) "
2357  << std::fabs(sum - 1.0) << " for node number " << n
2358  << " at value " << i << std::endl;
2359  }
2360  }
2361  }
2362  }
2363 #endif
2364  }
oomph::TreeBasedRefineableMeshBase::complete_hanging_nodes_recursively
void complete_hanging_nodes_recursively(Node *&nod_pt, Vector< Node * > &master_nodes, Vector< double > &hang_weights, const int &ival)
Auxiliary routine for recursive hanging node completion.
Definition: refineable_mesh.cc:2214
k
char char char int int * k
Definition: level2_impl.h:374
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117

References complete_hanging_nodes_recursively(), boost::multiprecision::fabs(), oomph::Node::hanging_pt(), i, oomph::Node::is_hanging(), k, oomph::HangInfo::master_weight(), n, oomph::HangInfo::nmaster(), oomph::Mesh::nnode(), oomph::Mesh::node_pt(), oomph::oomph_info, oomph::Node::set_hanging_pt(), and oomph::HangInfo::set_master_node_pt().

Referenced by adapt_mesh(), and p_adapt_mesh().

◆ complete_hanging_nodes_recursively()

void oomph::TreeBasedRefineableMeshBase::complete_hanging_nodes_recursively ( Node *&  nod_pt,
Vector< Node * > &  master_nodes,
Vector< double > &  hang_weights,
const int &  i 
)
protected

Auxiliary routine for recursive hanging node completion.

Given a node, return a vector of pointers to master nodes and a vector of the associated weights. This is done recursively, so if a node is not hanging, the node is regarded as its own master node which has weight 1.0.

2219  {
2220  // Is the node hanging in the variable i
2221  if (nod_pt->is_hanging(i))
2222  {
2223  // Loop over all master nodes
2224  HangInfo* const hang_pt = nod_pt->hanging_pt(i);
2225  unsigned nmaster = hang_pt->nmaster();
2226 
2227  for (unsigned m = 0; m < nmaster; m++)
2228  {
2229  // Get the master node
2230  Node* master_nod_pt = hang_pt->master_node_pt(m);
2231 
2232  // Keep in memory the size of the list before adding the nodes this
2233  // master node depends on. This is required so that the recursion is
2234  // only performed on these particular master nodes. A master node
2235  // could contain contributions from two separate pseudo-masters.
2236  // These contributions must be summed, not multiplied.
2237  int first_new_node = master_nodes.size();
2238 
2239  // Now check which master nodes this master node depends on
2240  complete_hanging_nodes_recursively(
2241  master_nod_pt, master_nodes, hang_weights, i);
2242 
2243  // Multiply old weight by new weight for all the nodes this master
2244  // node depends on
2245  unsigned n_new_master_node = master_nodes.size();
2246 
2247  double mtr_weight = hang_pt->master_weight(m);
2248 
2249  for (unsigned k = first_new_node; k < n_new_master_node; k++)
2250  {
2251  hang_weights[k] = mtr_weight * hang_weights[k];
2252  }
2253  }
2254  }
2255  else
2256  // Node isn't hanging so it enters itself with the full weight
2257  {
2258  master_nodes.push_back(nod_pt);
2259  hang_weights.push_back(1.0);
2260  }
2261  }
m
int * m
Definition: level2_cplx_impl.h:294

References oomph::Node::hanging_pt(), i, oomph::Node::is_hanging(), k, m, oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), and oomph::HangInfo::nmaster().

Referenced by complete_hanging_nodes().

◆ doc_adaptivity_targets()

void oomph::TreeBasedRefineableMeshBase::doc_adaptivity_targets ( std::ostream &  outfile)
inlinevirtual

Doc the targets for mesh adaptation.

Reimplemented from oomph::RefineableMeshBase.

468  {
469  outfile << std::endl;
470  outfile << "Targets for mesh adaptation: " << std::endl;
471  outfile << "---------------------------- " << std::endl;
472  outfile << "Target for max. error: " << Max_permitted_error << std::endl;
473  outfile << "Target for min. error: " << Min_permitted_error << std::endl;
474  outfile << "Min. refinement level: " << Min_refinement_level << std::endl;
475  outfile << "Max. refinement level: " << Max_refinement_level << std::endl;
476  outfile << "Min. p-refinement level: " << Min_p_refinement_level
477  << std::endl;
478  outfile << "Max. p-refinement level: " << Max_p_refinement_level
479  << std::endl;
480  outfile << "Don't unrefine if less than " << Max_keep_unrefined
481  << " elements need unrefinement." << std::endl;
482  outfile << std::endl;
483  }
oomph::RefineableMeshBase::Max_keep_unrefined
unsigned Max_keep_unrefined
Definition: refineable_mesh.h:359
oomph::RefineableMeshBase::Min_permitted_error
double Min_permitted_error
Min. error (i.e. (try to) merge elements if their error is smaller)
Definition: refineable_mesh.h:329
oomph::RefineableMeshBase::Max_permitted_error
double Max_permitted_error
Max. error (i.e. split elements if their error is larger)
Definition: refineable_mesh.h:326

References oomph::RefineableMeshBase::Max_keep_unrefined, Max_p_refinement_level, oomph::RefineableMeshBase::Max_permitted_error, Max_refinement_level, Min_p_refinement_level, oomph::RefineableMeshBase::Min_permitted_error, and Min_refinement_level.

Referenced by FSIRingProblem::dynamic_run().

◆ dump_refinement()

void oomph::TreeBasedRefineableMeshBase::dump_refinement ( std::ostream &  outfile)
virtual

Dump refinement pattern to allow for rebuild.

Dump refinement pattern to allow for rebuild

219  {
220  // Assign storage for refinement pattern
221  Vector<Vector<unsigned>> to_be_refined;
222 
223  // Get refinement pattern of reference mesh:
224  get_refinement_pattern(to_be_refined);
225 
226  // Dump max refinement level:
227  unsigned max_level = to_be_refined.size();
228  outfile << max_level << " # max. refinement level " << std::endl;
229 
230  // Doc the numbers of the elements that need to be refined at this level
231  for (unsigned l = 0; l < max_level; l++)
232  {
233  // Loop over elements that need to be refined at this level
234  unsigned n_to_be_refined = to_be_refined[l].size();
235  outfile << n_to_be_refined << " # number of elements to be refined. "
236  << "What follows are the numbers of the elements. " << std::endl;
237 
238  // Select relevant elements to be refined
239  for (unsigned i = 0; i < n_to_be_refined; i++)
240  {
241  outfile << to_be_refined[l][i] << std::endl;
242  }
243  }
244  }
oomph::TreeBasedRefineableMeshBase::get_refinement_pattern
virtual void get_refinement_pattern(Vector< Vector< unsigned >> &to_be_refined)
Definition: refineable_mesh.cc:50

References get_refinement_pattern(), and i.

◆ forest_pt()

TreeForest* oomph::TreeBasedRefineableMeshBase::forest_pt ( )
inline

Return pointer to the Forest represenation of the mesh.

461  {
462  return Forest_pt;
463  }

References Forest_pt.

Referenced by get_elements_at_refinement_level(), get_refinement_pattern(), p_adapt_mesh(), and refine_as_in_reference_mesh().

◆ get_elements_at_refinement_level()

void oomph::TreeBasedRefineableMeshBase::get_elements_at_refinement_level ( unsigned &  refinement_level,
Vector< RefineableElement * > &  level_elements 
)
virtual

Extract the elements at a particular refinement level in the refinement pattern - used in Mesh::redistribute or whatever it's going to be called (RefineableMeshBase::reduce_halo_layers or something)

Extract the elements at a particular refinement level in the refinement pattern (used in Mesh::redistribute or whatever it's going to be called (RefineableMeshBase::reduce_halo_layers or something)

115  {
116  // Extract *all* elements from current (fully refined) mesh.
117  Vector<Tree*> all_tree_nodes_pt;
118  forest_pt()->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
119 
120  // Add the element to the vector if its level matches refinement_level
121  unsigned nnodes = all_tree_nodes_pt.size();
122  for (unsigned e = 0; e < nnodes; e++)
123  {
124  unsigned level = all_tree_nodes_pt[e]->level();
125  if (level == refinement_level)
126  {
127  level_elements.push_back(
128  dynamic_cast<RefineableElement*>(all_tree_nodes_pt[e]->object_pt()));
129  }
130  }
131  }
oomph::TreeBasedRefineableMeshBase::forest_pt
TreeForest * forest_pt()
Return pointer to the Forest represenation of the mesh.
Definition: refineable_mesh.h:460

References e(), forest_pt(), and oomph::TreeForest::stick_all_tree_nodes_into_vector().

◆ get_refinement_levels()

void oomph::TreeBasedRefineableMeshBase::get_refinement_levels ( unsigned &  min_refinement_level,
unsigned &  max_refinement_level 
)
virtual

Get max/min refinement levels in mesh.

Get max/min refinement level.

832  {
833  // Initialise
834  min_refinement_level = UINT_MAX;
835  max_refinement_level = 0;
836 
837  // Loop over all elements
838  unsigned long n_element = this->nelement();
839  if (n_element == 0)
840  {
841  min_refinement_level = 0;
842  max_refinement_level = 0;
843  }
844  else
845  {
846  for (unsigned long e = 0; e < n_element; e++)
847  {
848  // Get the refinement level of the element
849  unsigned level = dynamic_cast<RefineableElement*>(this->element_pt(e))
850  ->refinement_level();
851 
852  if (level > max_refinement_level) max_refinement_level = level;
853  if (level < min_refinement_level) min_refinement_level = level;
854  }
855  }
856  }

References e(), oomph::Mesh::element_pt(), max_refinement_level(), min_refinement_level(), and oomph::Mesh::nelement().

Referenced by oomph::Problem::read(), refine_as_in_reference_mesh(), refine_base_mesh(), oomph::RefineableBrickMesh< ELEMENT >::setup_octree_forest(), and oomph::RefineableQuadMesh< ELEMENT >::setup_quadtree_forest().

◆ get_refinement_pattern()

void oomph::TreeBasedRefineableMeshBase::get_refinement_pattern ( Vector< Vector< unsigned >> &  to_be_refined)
virtual

Extract refinement pattern: Consider the hypothetical mesh obtained by truncating the refinement of the current mesh to a given level (where level=0 is the un-refined base mesh). To advance to the next refinement level, we need to refine (split) the to_be_refined[level].size() elements identified by the element numbers contained in vector to_be_refined[level][...]

Get refinement pattern of mesh: Consider the hypothetical mesh obtained by truncating the refinement of the current mesh to a given level (where level=0 is the un-refined base mesh). To advance to the next refinement level, we need to refine (split) the to_be_refined[level].size() elements identified by the element numbers contained in vector to_be_refined[level][...]

52  {
53  // Extract *all* elements from current (fully refined) mesh.
54  Vector<Tree*> all_tree_nodes_pt;
55  forest_pt()->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
56 
57  // Find out maximum refinement level
58  unsigned max_level = 0;
59  unsigned nnodes = all_tree_nodes_pt.size();
60  for (unsigned e = 0; e < nnodes; e++)
61  {
62  unsigned level = all_tree_nodes_pt[e]->level();
63  if (level > max_level) max_level = level;
64  }
65 
66  // Assign storage for refinement pattern
67  to_be_refined.clear();
68  to_be_refined.resize(max_level);
69  Vector<unsigned> el_count(max_level);
70 
71  // Initialise count of elements that exist in mesh when refinement
72  // has proceeded to this level
73  for (unsigned l = 0; l < max_level; l++)
74  {
75  el_count[l] = 0;
76  }
77 
78  // Loop over all levels and extract all elements that exist
79  // in reference mesh when refinement has proceeded to this level
80  for (unsigned l = 0; l < max_level; l++)
81  {
82  // Loop over all elements (tree nodes)
83  for (unsigned e = 0; e < nnodes; e++)
84  {
85  // What level does this element exist on?
86  unsigned level = all_tree_nodes_pt[e]->level();
87 
88  // Element is part of the mesh at this refinement level
89  // if it exists at this level OR if it exists at a lower level
90  // and is a leaf
91  if ((level == l) || ((level < l) && (all_tree_nodes_pt[e]->is_leaf())))
92  {
93  // If element exsts at this level and is not a leaf it will
94  // be refined when we move to the next level:
95  if ((level == l) && (!all_tree_nodes_pt[e]->is_leaf()))
96  {
97  // Add element number (in mesh at current refinement level)
98  // to the list of elements that need to be refined
99  to_be_refined[l].push_back(el_count[l]);
100  }
101  // Element exists in this mesh: Add to counter
102  el_count[l]++;
103  }
104  }
105  }
106  }

References e(), forest_pt(), and oomph::TreeForest::stick_all_tree_nodes_into_vector().

Referenced by dump_refinement(), refine_base_mesh_as_in_reference_mesh(), and refine_base_mesh_as_in_reference_mesh_minus_one().

◆ max_p_refinement_level()

unsigned& oomph::TreeBasedRefineableMeshBase::max_p_refinement_level ( )
inline

Access fct for max. permissible p-refinement level (relative to base mesh)

500  {
501  return Max_p_refinement_level;
502  }

References Max_p_refinement_level.

Referenced by p_adapt().

◆ max_refinement_level()

unsigned& oomph::TreeBasedRefineableMeshBase::max_refinement_level ( )
inline

Access fct for max. permissible refinement level (relative to base mesh)

487  {
488  return Max_refinement_level;
489  }

References Max_refinement_level.

Referenced by adapt(), FSIRingProblem::dynamic_run(), and get_refinement_levels().

◆ min_p_refinement_level()

unsigned& oomph::TreeBasedRefineableMeshBase::min_p_refinement_level ( )
inline

Access fct for min. permissible p-refinement level (relative to base mesh)

507  {
508  return Min_p_refinement_level;
509  }

References Min_p_refinement_level.

◆ min_refinement_level()

unsigned& oomph::TreeBasedRefineableMeshBase::min_refinement_level ( )
inline

Access fct for min. permissible refinement level (relative to base mesh)

493  {
494  return Min_refinement_level;
495  }

References Min_refinement_level.

Referenced by adapt(), FSIRingProblem::dynamic_run(), and get_refinement_levels().

◆ operator=()

void oomph::TreeBasedRefineableMeshBase::operator= ( const TreeBasedRefineableMeshBase &  )
delete

Broken assignment operator.

◆ p_adapt()

void oomph::TreeBasedRefineableMeshBase::p_adapt ( const Vector< double > &  elemental_error)
virtual

p-adapt mesh: Refine elements whose error is lager than err_max and (try to) unrefine those whose error is smaller than err_min

///////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////// Do adaptive p-refinement for mesh.

  • Pass Vector of error estimates for all elements.
  • p-refine those whose errors exceeds the threshold
  • p-unrefine those whose errors is less than threshold.

/ ... otherwise mark it as having been over-ruled

Reimplemented from oomph::RefineableMeshBase.

3953  {
3954  // Set the refinement tolerance to be the max permissible error
3955  double refine_tol = this->max_permitted_error();
3956 
3957  // Set the unrefinement tolerance to be the min permissible error
3958  double unrefine_tol = this->min_permitted_error();
3959 
3960  // Setup doc info
3961  DocInfo local_doc_info;
3962  if (doc_info_pt() == 0)
3963  {
3964  local_doc_info.disable_doc();
3965  }
3966  else
3967  {
3968  local_doc_info = this->doc_info();
3969  }
3970 
3971 
3972  // Check that the errors make sense
3973  if (refine_tol <= unrefine_tol)
3974  {
3975  std::ostringstream error_stream;
3976  error_stream << "Refinement tolerance <= Unrefinement tolerance"
3977  << refine_tol << " " << unrefine_tol << std::endl
3978  << "doesn't make sense and will almost certainly crash"
3979  << std::endl
3980  << "this beautiful code!" << std::endl;
3981 
3982  throw OomphLibError(
3983  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3984  }
3985 
3986 
3987  // Select elements for refinement and unrefinement
3988  //==============================================
3989  // Reset counter for number of elements that would like to be
3990  // refined further but can't
3991  this->nrefinement_overruled() = 0;
3992 
3993  unsigned n_refine = 0;
3994  unsigned n_unrefine = 0;
3995  // Loop over all elements and mark them according to the error criterion
3996  unsigned long Nelement = this->nelement();
3997  for (unsigned long e = 0; e < Nelement; e++)
3998  {
3999  //(Cast) pointer to the element
4000  PRefineableElement* el_pt =
4001  dynamic_cast<PRefineableElement*>(this->element_pt(e));
4002 
4003  // Check that we can p-refine the element
4004  if (el_pt != 0)
4005  {
4006  // Initially element is not to be refined
4007  el_pt->deselect_for_p_refinement();
4008  el_pt->deselect_for_p_unrefinement();
4009 
4010  // If the element error exceeds the threshold ...
4011  if (elemental_error[e] > refine_tol)
4012  {
4013  // ... and its refinement level is less than the maximum desired level
4014  // mark is to be refined
4015  if ((el_pt->p_refinement_is_enabled()) &&
4016  (el_pt->p_order() < this->max_p_refinement_level()))
4017  {
4018  el_pt->select_for_p_refinement();
4019  n_refine++;
4020  }
4021  // ... otherwise mark it as having been over-ruled
4022  else
4023  {
4024  this->nrefinement_overruled() += 1;
4025  }
4026  }
4027  if (elemental_error[e] < unrefine_tol)
4028  {
4029  // ... and its refinement level is more than the minimum desired level
4030  // mark is to be refined
4031  //(Also check that we don't unrefine past the initial refinement
4032  // level)
4033  if ((el_pt->p_refinement_is_enabled()) &&
4034  (el_pt->p_order() > this->min_p_refinement_level()) &&
4035  (el_pt->p_order() > el_pt->initial_p_order()))
4036  {
4037  el_pt->select_for_p_unrefinement();
4038  n_unrefine++;
4039  }
4040  // Don't mark as overruled - it's misleading
4042  // else
4043  // {
4044  // this->nrefinement_overruled()+=1;
4045  // }
4046  }
4047  } // End of check for p-refineability of element
4048  else
4049  {
4050  oomph_info << "p-refinement is not possible for these elements"
4051  << std::endl;
4052  // Don't try to p-refine any more elements
4053  break;
4054  }
4055  }
4056 
4057  oomph_info << " \n Number of elements to be refined: " << n_refine
4058  << std::endl;
4059  oomph_info << " \n Number of elements whose refinement was overruled: "
4060  << this->nrefinement_overruled() << std::endl;
4061 
4062  oomph_info << " \n Number of elements to be unrefined : " << n_unrefine
4063  << std::endl
4064  << std::endl;
4065 
4066 
4067  // Now do the actual mesh adaptation
4068  //---------------------------------
4069 
4070  // Check whether its worth our while
4071  // Either some elements want to be refined,
4072  // or the number that want to be unrefined are greater than the
4073  // specified tolerance
4074 
4075  // In a parallel job, it is possible that one process may not have
4076  // any elements to refine, BUT a neighbouring process may refine an
4077  // element which changes the hanging status of a node that is on
4078  // both processes (i.e. a halo(ed) node). To get around this issue,
4079  // ALL processes need to call adapt_mesh if ANY refinement is to
4080  // take place anywhere.
4081 
4082  unsigned total_n_refine = 0;
4083 #ifdef OOMPH_HAS_MPI
4084  // Sum n_refine across all processors
4085  if (this->is_mesh_distributed())
4086  {
4087  MPI_Allreduce(
4088  &n_refine, &total_n_refine, 1, MPI_INT, MPI_SUM, Comm_pt->mpi_comm());
4089  }
4090  else
4091  {
4092  total_n_refine = n_refine;
4093  }
4094 #else
4095  total_n_refine = n_refine;
4096 #endif
4097 
4098  // There may be some issues with unrefinement too, but I have not
4099  // been able to come up with an example (either in my head or in a
4100  // particular problem) where anything has arisen. I can see that
4101  // there may be an issue if n_unrefine differs across processes so
4102  // that (total_n_unrefine > max_keep_unrefined()) on some but not
4103  // all processes. I haven't seen any examples of this yet so the
4104  // following code may or may not work! (Andy, 06/03/08)
4105 
4106  unsigned total_n_unrefine = 0;
4107 #ifdef OOMPH_HAS_MPI
4108  // Sum n_unrefine across all processors
4109  if (this->is_mesh_distributed())
4110  {
4111  MPI_Allreduce(&n_unrefine,
4112  &total_n_unrefine,
4113  1,
4114  MPI_INT,
4115  MPI_SUM,
4116  Comm_pt->mpi_comm());
4117  }
4118  else
4119  {
4120  total_n_unrefine = n_unrefine;
4121  }
4122 #else
4123  total_n_unrefine = n_unrefine;
4124 #endif
4125 
4126  oomph_info << "---> " << total_n_refine << " elements to be refined, and "
4127  << total_n_unrefine << " to be unrefined, in total."
4128  << std::endl;
4129 
4130  if ((total_n_refine > 0) || (total_n_unrefine > this->max_keep_unrefined()))
4131  {
4132 #ifdef PARANOID
4133 #ifdef OOMPH_HAS_MPI
4134 
4135  // Sanity check: Each processor checks if the enforced unrefinement of
4136  // its haloed element is matched by enforced unrefinement of the
4137  // corresponding halo elements on the other processors.
4138  if (this->is_mesh_distributed())
4139  {
4140  // Store number of processors and current process
4141  MPI_Status status;
4142  int n_proc = Comm_pt->nproc();
4143  int my_rank = Comm_pt->my_rank();
4144 
4145  // Loop over all other domains/processors
4146  for (int d = 0; d < n_proc; d++)
4147  {
4148  // Don't talk to yourself
4149  if (d != my_rank)
4150  {
4151  {
4152  // Get the vector of halo elements whose non-halo counterpart
4153  // are on processor d
4154  Vector<GeneralisedElement*> halo_elem_pt(
4155  this->halo_element_pt(d));
4156 
4157  // Create vector containing (0)1 to indicate that
4158  // halo element is (not) to be unrefined
4159  unsigned nhalo = halo_elem_pt.size();
4160  Vector<int> halo_to_be_unrefined(nhalo, 0);
4161  for (unsigned e = 0; e < nhalo; e++)
4162  {
4163  if (dynamic_cast<PRefineableElement*>(halo_elem_pt[e])
4164  ->to_be_p_unrefined())
4165  {
4166  halo_to_be_unrefined[e] = 1;
4167  }
4168  }
4169 
4170  // Trap the case when there are no halo elements
4171  // so that we don't get a segfault in the MPI send
4172  if (nhalo > 0)
4173  {
4174  // Send it across
4175  MPI_Send(&halo_to_be_unrefined[0],
4176  nhalo,
4177  MPI_INT,
4178  d,
4179  0,
4180  Comm_pt->mpi_comm());
4181  }
4182  }
4183 
4184  {
4185  // Get the vector of haloed elements on current processor
4186  Vector<GeneralisedElement*> haloed_elem_pt(
4187  this->haloed_element_pt(d));
4188 
4189  // Ask processor d to send vector containing (0)1 for
4190  // halo element with current processor to be (not)unrefined
4191  unsigned nhaloed = haloed_elem_pt.size();
4192  Vector<int> halo_to_be_unrefined(nhaloed);
4193  // Trap to catch the case that there are no haloed elements
4194  if (nhaloed > 0)
4195  {
4196  MPI_Recv(&halo_to_be_unrefined[0],
4197  nhaloed,
4198  MPI_INT,
4199  d,
4200  0,
4201  Comm_pt->mpi_comm(),
4202  &status);
4203  }
4204 
4205  // Check it
4206  for (unsigned e = 0; e < nhaloed; e++)
4207  {
4208  if (((halo_to_be_unrefined[e] == 0) &&
4209  (dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])
4210  ->to_be_p_unrefined())) ||
4211  ((halo_to_be_unrefined[e] == 1) &&
4212  (!dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])
4213  ->to_be_p_unrefined())))
4214  {
4215  std::ostringstream error_message;
4216  error_message
4217  << "Error in refinement: \n"
4218  << "Haloed element: " << e << " on proc " << my_rank
4219  << " \n"
4220  << "wants to be unrefined whereas its halo counterpart on\n"
4221  << "proc " << d << " doesn't (or vice versa)...\n"
4222  << "This is most likely because the error estimator\n"
4223  << "has not assigned the same errors to halo and haloed\n"
4224  << "elements -- it ought to!\n";
4225  throw OomphLibError(error_message.str(),
4226  OOMPH_CURRENT_FUNCTION,
4227  OOMPH_EXCEPTION_LOCATION);
4228  }
4229  }
4230  }
4231  }
4232  }
4233 
4234 
4235  // Loop over all other domains/processors
4236  for (int d = 0; d < n_proc; d++)
4237  {
4238  // Don't talk to yourself
4239  if (d != my_rank)
4240  {
4241  {
4242  // Get the vector of halo elements whose non-halo counterpart
4243  // are on processor d
4244  Vector<GeneralisedElement*> halo_elem_pt(
4245  this->halo_element_pt(d));
4246 
4247  // Create vector containing (0)1 to indicate that
4248  // halo element is (not) to be refined
4249  unsigned nhalo = halo_elem_pt.size();
4250  Vector<int> halo_to_be_refined(nhalo, 0);
4251  for (unsigned e = 0; e < nhalo; e++)
4252  {
4253  if (dynamic_cast<PRefineableElement*>(halo_elem_pt[e])
4254  ->to_be_p_refined())
4255  {
4256  halo_to_be_refined[e] = 1;
4257  }
4258  }
4259 
4260  // Send it across
4261  if (nhalo > 0)
4262  {
4263  MPI_Send(&halo_to_be_refined[0],
4264  nhalo,
4265  MPI_INT,
4266  d,
4267  0,
4268  Comm_pt->mpi_comm());
4269  }
4270  }
4271 
4272  {
4273  // Get the vector of haloed elements on current processor
4274  Vector<GeneralisedElement*> haloed_elem_pt(
4275  this->haloed_element_pt(d));
4276 
4277  // Ask processor d to send vector containing (0)1 for
4278  // halo element with current processor to be (not)refined
4279  unsigned nhaloed = haloed_elem_pt.size();
4280  Vector<int> halo_to_be_refined(nhaloed);
4281  if (nhaloed > 0)
4282  {
4283  MPI_Recv(&halo_to_be_refined[0],
4284  nhaloed,
4285  MPI_INT,
4286  d,
4287  0,
4288  Comm_pt->mpi_comm(),
4289  &status);
4290  }
4291 
4292  // Check it
4293  for (unsigned e = 0; e < nhaloed; e++)
4294  {
4295  if (((halo_to_be_refined[e] == 0) &&
4296  (dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])
4297  ->to_be_p_refined())) ||
4298  ((halo_to_be_refined[e] == 1) &&
4299  (!dynamic_cast<PRefineableElement*>(haloed_elem_pt[e])
4300  ->to_be_p_refined())))
4301  {
4302  std::ostringstream error_message;
4303  error_message
4304  << "Error in refinement: \n"
4305  << "Haloed element: " << e << " on proc " << my_rank
4306  << " \n"
4307  << "wants to be refined whereas its halo counterpart on\n"
4308  << "proc " << d << " doesn't (or vice versa)...\n"
4309  << "This is most likely because the error estimator\n"
4310  << "has not assigned the same errors to halo and haloed\n"
4311  << "elements -- it ought to!\n";
4312  throw OomphLibError(error_message.str(),
4313  OOMPH_CURRENT_FUNCTION,
4314  OOMPH_EXCEPTION_LOCATION);
4315  }
4316  }
4317  }
4318  }
4319  }
4320  }
4321 #endif
4322 #endif
4323 
4324  // Perform the actual adaptation
4325  p_adapt_mesh(local_doc_info);
4326 
4327  // The number of refineable elements is still local to each process
4328  this->Nunrefined = n_unrefine;
4329  this->Nrefined = n_refine;
4330  }
4331  // If not worthwhile, say so but still reorder nodes and kill external
4332  // storage for consistency in parallel computations
4333  else
4334  {
4335 #ifdef OOMPH_HAS_MPI
4336  // Delete any external element storage - any interaction will still
4337  // be set up on the fly again, so we need to get rid of old information.
4338  // This particularly causes problems in multi-domain examples where
4339  // we decide not to refine one of the meshes
4340  this->delete_all_external_storage();
4341 #endif
4342 
4343  // Reorder the nodes within the mesh's node vector
4344  // to establish a standard ordering regardless of the sequence
4345  // of mesh refinements -- this is required to allow dump/restart
4346  // on refined meshes
4347  this->reorder_nodes();
4348 
4349 #ifdef OOMPH_HAS_MPI
4350 
4351  // Now (re-)classify halo and haloed nodes and synchronise hanging
4352  // nodes
4353  // This is required in cases where delete_all_external_storage()
4354  // made dependent nodes to external halo nodes nonhanging.
4355  if (this->is_mesh_distributed())
4356  {
4357  DocInfo doc_info;
4358  doc_info.disable_doc();
4359  classify_halo_and_haloed_nodes(doc_info, doc_info.is_doc_enabled());
4360  }
4361 
4362 #endif
4363 
4364  if (n_refine == 0)
4365  {
4366  oomph_info << "\n Not enough benefit in adapting mesh. " << std::endl
4367  << std::endl;
4368  }
4369  this->Nunrefined = 0;
4370  this->Nrefined = 0;
4371  }
4372  }
oomph::TreeBasedRefineableMeshBase::max_p_refinement_level
unsigned & max_p_refinement_level()
Definition: refineable_mesh.h:499
oomph::TreeBasedRefineableMeshBase::p_adapt_mesh
void p_adapt_mesh()
Definition: refineable_mesh.h:533

References oomph::Mesh::delete_all_external_storage(), oomph::PRefineableElement::deselect_for_p_refinement(), oomph::PRefineableElement::deselect_for_p_unrefinement(), oomph::DocInfo::disable_doc(), oomph::RefineableMeshBase::doc_info(), oomph::RefineableMeshBase::doc_info_pt(), e(), oomph::Mesh::element_pt(), oomph::PRefineableElement::initial_p_order(), oomph::DocInfo::is_doc_enabled(), oomph::Mesh::is_mesh_distributed(), oomph::RefineableMeshBase::max_keep_unrefined(), max_p_refinement_level(), oomph::RefineableMeshBase::max_permitted_error(), oomph::RefineableMeshBase::min_permitted_error(), oomph::Mesh::nelement(), oomph::RefineableMeshBase::Nrefined, oomph::RefineableMeshBase::nrefinement_overruled(), oomph::RefineableMeshBase::Nunrefined, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, p_adapt_mesh(), oomph::PRefineableElement::p_order(), oomph::PRefineableElement::p_refinement_is_enabled(), oomph::Mesh::reorder_nodes(), oomph::PRefineableElement::select_for_p_refinement(), and oomph::PRefineableElement::select_for_p_unrefinement().

◆ p_adapt_mesh() [1/2]

void oomph::TreeBasedRefineableMeshBase::p_adapt_mesh ( )
inline

Perform the actual tree-based mesh p-adaptation. A simple wrapper to call the function without documentation.

534  {
535  // Create a dummy doc_info object
536  DocInfo doc_info;
537  doc_info.directory() = "";
538  doc_info.disable_doc();
539  // Call the other adapt mesh
540  p_adapt_mesh(doc_info);
541  }

References oomph::DocInfo::directory(), oomph::DocInfo::disable_doc(), and oomph::RefineableMeshBase::doc_info().

Referenced by p_adapt(), p_refine_selected_elements(), p_refine_uniformly(), and p_unrefine_uniformly().

◆ p_adapt_mesh() [2/2]

void oomph::TreeBasedRefineableMeshBase::p_adapt_mesh ( DocInfo &  doc_info)

Perform the actual tree-based mesh p-adaptation, documenting the progress in the directory specified in DocInfo object.

p-adapt mesh, which exists in two representations, namely as:

  • a FE mesh
  • a forest of Oc or QuadTrees

p-refinement/derefinement process is documented (in tecplot-able form) if requested.

Procedure:

  • Loop over all elements and do the p-refinement/unrefinement for those who want to be refined. Note: p-refinement builds fully- functional elements.
  • For all nodes that were hanging on the previous mesh (and are still marked as such), fill in their nodal values (consistent with the current hanging node scheme) to make sure they are fully functional, should they have become non-hanging during the mesh-adaptation. Then mark the nodes as non-hanging.
  • Delete any nodes that have become obsolete.
  • Mark up hanging nodes and setup hanging node scheme (incl. recursive cleanup for hanging nodes that depend on other hanging nodes).
  • run a quick self-test on the neighbour finding scheme and check the integrity of the elements (if PARANOID)
  • doc hanging node status, boundary conditions, neighbour scheme if requested.

After adaptation, all nodes (whether new or old) have up-to-date current and previous values.

If refinement process is being documented, the following information is documented:

  • The files
    • "neighbours.dat"
    • "all_nodes.dat"
    • "new_nodes.dat"
    • "hang_nodes_*.dat" where the * denotes a direction (n,s,e,w) in 2D or (r,l,u,d,f,b) in 3D

can be viewed with

  • QHangingNodes.mcr
  • The file
    • "hangnodes_withmasters.dat"

can be viewed with

  • QHangingNodesWithMasters.mcr

to check the hanging node status.

  • The neighbour status of the elements is documented in
    • "neighbours.dat"

and can be viewed with

  • QuadTreeNeighbours.mcr

Update the boundary element info – this can be a costly procedure and for this reason the mesh writer might have decided not to set up this scheme. If so, we won't change this and suppress its creation...

/BENFLAG: Check that all the nodes belong to their update elements

4434  {
4435 #ifdef OOMPH_HAS_MPI
4436  // Delete any external element storage before performing the adaptation
4437  // (in particular, external halo nodes that are on mesh boundaries)
4438  this->delete_all_external_storage();
4439 #endif
4440 
4441  // Only perform the adapt step if the mesh has any elements. This is
4442  // relevant in a distributed problem with multiple meshes, where a
4443  // particular process may not have any elements on a particular submesh.
4444  if (this->nelement() > 0)
4445  {
4446  double t_start = 0.0;
4447  if (Global_timings::Doc_comprehensive_timings)
4448  {
4449  t_start = TimingHelpers::timer();
4450  }
4451 
4452  // Do refinement/unrefinement if required
4453  this->p_refine_elements_if_required();
4454 
4455  if (Global_timings::Doc_comprehensive_timings)
4456  {
4457  double t_end = TimingHelpers::timer();
4458  oomph_info << "Time for p-refinement/unrefinement: " << t_end - t_start
4459  << std::endl;
4460  t_start = TimingHelpers::timer();
4461  }
4462 
4463  // Loop over all nodes in mesh and free the dofs of those that were
4464  //-----------------------------------------------------------------
4465  // pinned only because they were hanging nodes. Also update their
4466  //-----------------------------------------------------------------
4467  // nodal values so that they contain data that is consistent
4468  //----------------------------------------------------------
4469  // with the hanging node representation
4470  //-------------------------------------
4471  // (Even if the nodal data isn't actually accessed because the node
4472  // is still hanging -- we don't know this yet, and this step makes
4473  // sure that all nodes are fully functional and up-to-date, should
4474  // they become non-hanging below).
4475  //
4476  //
4477  // However, if we have a fixed mesh and hanging nodes on the boundary
4478  // become non-hanging they will not necessarily respect the curvilinear
4479  // boundaries. This can only happen in 3D of course because it is not
4480  // possible to have hanging nodes on boundaries in 2D.
4481  // The solution is to store those nodes on the boundaries that are
4482  // currently hanging and then check to see whether they have changed
4483  // status at the end of the refinement procedure.
4484  // If it has changed, then we need to adjust their positions.
4485  const unsigned n_boundary = this->nboundary();
4486  const unsigned mesh_dim = this->finite_element_pt(0)->dim();
4487  Vector<std::set<Node*>> hanging_nodes_on_boundary_pt(n_boundary);
4488 
4489  unsigned long n_node = this->nnode();
4490  for (unsigned long n = 0; n < n_node; n++)
4491  {
4492  // Get the pointer to the node
4493  Node* nod_pt = this->node_pt(n);
4494 
4495  // Get the number of values in the node
4496  unsigned n_value = nod_pt->nvalue();
4497 
4498  // We need to find if any of the values are hanging
4499  bool is_hanging = nod_pt->is_hanging();
4500  // Loop over the values and find out whether any are hanging
4501  for (unsigned n = 0; n < n_value; n++)
4502  {
4503  is_hanging |= nod_pt->is_hanging(n);
4504  }
4505 
4506  // If the node is hanging then ...
4507  if (is_hanging)
4508  {
4509  // Unless they are turned into hanging nodes again below
4510  // (this might or might not happen), fill in all the necessary
4511  // data to make them 'proper' nodes again.
4512 
4513  // Reconstruct the nodal values/position from the node's
4514  // hanging node representation
4515  unsigned nt = nod_pt->ntstorage();
4516  Vector<double> values(n_value);
4517  unsigned n_dim = nod_pt->ndim();
4518  Vector<double> position(n_dim);
4519  // Loop over all history values
4520  for (unsigned t = 0; t < nt; t++)
4521  {
4522  nod_pt->value(t, values);
4523  for (unsigned i = 0; i < n_value; i++)
4524  {
4525  nod_pt->set_value(t, i, values[i]);
4526  }
4527  nod_pt->position(t, position);
4528  for (unsigned i = 0; i < n_dim; i++)
4529  {
4530  nod_pt->x(t, i) = position[i];
4531  }
4532  }
4533 
4534  // If it's an algebraic node: Update its previous nodal positions too
4535  AlgebraicNode* alg_node_pt = dynamic_cast<AlgebraicNode*>(nod_pt);
4536  if (alg_node_pt != 0)
4537  {
4538  bool update_all_time_levels = true;
4539  alg_node_pt->node_update(update_all_time_levels);
4540  }
4541 
4542 
4543  // If it's a Solid node, update Lagrangian coordinates
4544  // from its hanging node representation
4545  SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
4546  if (solid_node_pt != 0)
4547  {
4548  unsigned n_lagrangian = solid_node_pt->nlagrangian();
4549  for (unsigned i = 0; i < n_lagrangian; i++)
4550  {
4551  solid_node_pt->xi(i) = solid_node_pt->lagrangian_position(i);
4552  }
4553  }
4554 
4555  // Now store geometrically hanging nodes on boundaries that
4556  // may need updating after refinement.
4557  // There will only be a problem if we have 3 spatial dimensions
4558  if ((mesh_dim > 2) && (nod_pt->is_hanging()))
4559  {
4560  // If the node is on a boundary then add a pointer to the node
4561  // to our lookup scheme
4562  if (nod_pt->is_on_boundary())
4563  {
4564  // Storage for the boundaries on which the Node is located
4565  std::set<unsigned>* boundaries_pt;
4566  nod_pt->get_boundaries_pt(boundaries_pt);
4567  if (boundaries_pt != 0)
4568  {
4569  // Loop over the boundaries and add a pointer to the node
4570  // to the appropriate storage scheme
4571  for (std::set<unsigned>::iterator it = boundaries_pt->begin();
4572  it != boundaries_pt->end();
4573  ++it)
4574  {
4575  hanging_nodes_on_boundary_pt[*it].insert(nod_pt);
4576  }
4577  }
4578  }
4579  }
4580 
4581  } // End of is_hanging
4582 
4583  // Initially mark all nodes as 'non-hanging' and `obsolete'
4584  nod_pt->set_nonhanging();
4585  nod_pt->set_obsolete();
4586  }
4587 
4588  double t_end = 0.0;
4589  if (Global_timings::Doc_comprehensive_timings)
4590  {
4591  t_end = TimingHelpers::timer();
4592  oomph_info << "Time for sorting out initial hanging status: "
4593  << t_end - t_start << std::endl;
4594  t_start = TimingHelpers::timer();
4595  }
4596 
4597  // Stick all elements into a vector
4598  Vector<Tree*> tree_nodes_pt;
4599  this->forest_pt()->stick_leaves_into_vector(tree_nodes_pt);
4600 
4601  // Copy the elements into the mesh Vector
4602  unsigned long num_tree_nodes = tree_nodes_pt.size();
4603  this->element_pt().resize(num_tree_nodes);
4604  for (unsigned long e = 0; e < num_tree_nodes; e++)
4605  {
4606  this->element_pt(e) = tree_nodes_pt[e]->object_pt();
4607 
4608  // Now loop over all nodes in element and mark them as non-obsolete
4609  FiniteElement* this_el_pt = this->finite_element_pt(e);
4610  unsigned n_node = this_el_pt->nnode(); // caching pre-loop
4611  for (unsigned n = 0; n < n_node; n++)
4612  {
4613  this_el_pt->node_pt(n)->set_non_obsolete();
4614  }
4615 
4616  // Mark up so that repeated refinements do not occur
4617  // (Required because refined element is the same element as the
4618  // original)
4619  PRefineableElement* cast_el_pt =
4620  dynamic_cast<PRefineableElement*>(this->element_pt(e));
4621  cast_el_pt->deselect_for_p_refinement();
4622  cast_el_pt->deselect_for_p_unrefinement();
4623  }
4624 
4625  // Cannot delete nodes that are still marked as obsolete
4626  // because they may still be required to assemble the hanging schemes
4627  //-------------------------------------------------------------------
4628 
4629  // Mark up hanging nodes
4630  //----------------------
4631 
4632  // Output streams for the hanging nodes
4633  Vector<std::ofstream*> hanging_output_files;
4634  // Setup the output files for hanging nodes, this must be called
4635  // precisely once for the forest. Note that the files will only
4636  // actually be opened if doc_info.Doc_flag is true
4637  this->forest_pt()->open_hanging_node_files(doc_info,
4638  hanging_output_files);
4639 
4640  for (unsigned long e = 0; e < num_tree_nodes; e++)
4641  {
4642  // Generic setup
4643  tree_nodes_pt[e]->object_pt()->setup_hanging_nodes(
4644  hanging_output_files);
4645  // Element specific setup
4646  tree_nodes_pt[e]->object_pt()->further_setup_hanging_nodes();
4647  }
4648 
4649  // Close the hanging node files and delete the memory allocated
4650  // for the streams
4651  this->forest_pt()->close_hanging_node_files(doc_info,
4652  hanging_output_files);
4653 
4654 
4655  if (Global_timings::Doc_comprehensive_timings)
4656  {
4657  t_end = TimingHelpers::timer();
4658  oomph_info << "Time for setup_hanging_nodes() and "
4659  "further_setup_hanging_nodes() for "
4660  << num_tree_nodes << " elements: " << t_end - t_start
4661  << std::endl;
4662  t_start = TimingHelpers::timer();
4663  }
4664 
4665  // Read out the number of continously interpolated values
4666  // from one of the elements (assuming it's the same in all elements)
4667  unsigned ncont_interpolated_values =
4668  tree_nodes_pt[0]->object_pt()->ncont_interpolated_values();
4669 
4670  // Complete the hanging nodes schemes by dealing with the
4671  // recursively hanging nodes
4672  this->complete_hanging_nodes(ncont_interpolated_values);
4673 
4674 
4675  if (Global_timings::Doc_comprehensive_timings)
4676  {
4677  t_end = TimingHelpers::timer();
4678  oomph_info << "Time for complete_hanging_nodes: " << t_end - t_start
4679  << std::endl;
4680  t_start = TimingHelpers::timer();
4681  }
4682 
4686  if (Lookup_for_elements_next_boundary_is_setup)
4687  {
4688  this->setup_boundary_element_info();
4689  }
4690 
4691  if (Global_timings::Doc_comprehensive_timings)
4692  {
4693  t_end = TimingHelpers::timer();
4694  oomph_info << "Time for boundary element info: " << t_end - t_start
4695  << std::endl;
4696  t_start = TimingHelpers::timer();
4697  }
4698 
4699  // BENFLAG: Reset all the node update elements.
4700  // This is necessary to prevent the following case: A node N is
4701  // shared between two elements, A and B. The update element for
4702  // the node is set to A, say. Element A is p-refined and now
4703  // nolonger has N as a node. However the node update element for N
4704  // is still A but the node doesn't exist in A.
4705  MacroElementNodeUpdateElementBase* first_macro_el_pt =
4706  dynamic_cast<MacroElementNodeUpdateElementBase*>(this->element_pt(0));
4707  if (first_macro_el_pt != 0)
4708  {
4709  // Now set the node update info elementwise
4710  for (unsigned e = 0; e < this->nelement(); e++)
4711  {
4712  // Cast to macro element
4713  MacroElementNodeUpdateElementBase* macro_el_pt =
4714  dynamic_cast<MacroElementNodeUpdateElementBase*>(
4715  this->element_pt(e));
4716  if (macro_el_pt != 0)
4717  {
4718  // Get vector of geometric objects from element (construct vector
4719  // via copy operation)
4720  Vector<GeomObject*> geom_object_pt(macro_el_pt->geom_object_pt());
4721 
4722  // (Re)set node update info for all the nodes in the element
4723  macro_el_pt->set_node_update_info(geom_object_pt);
4724  }
4725  }
4726  }
4727 
4728 #ifdef PARANOID
4729 
4730  // Doc/check the neighbours
4731  //-------------------------
4732  Vector<Tree*> all_tree_nodes_pt;
4733  this->forest_pt()->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
4734 
4735  // Check the neighbours
4736  this->forest_pt()->check_all_neighbours(doc_info);
4737 
4738  // Check the integrity of the elements
4739  // -----------------------------------
4740 
4741  // Loop over elements and get the elemental integrity
4742  double max_error = 0.0;
4743  for (unsigned long e = 0; e < num_tree_nodes; e++)
4744  {
4745  double max_el_error;
4746  tree_nodes_pt[e]->object_pt()->check_integrity(max_el_error);
4747  // If the elemental error is greater than our maximum error
4748  // reset the maximum
4749  if (max_el_error > max_error)
4750  {
4751  max_error = max_el_error;
4752  }
4753  }
4754 
4755  if (max_error > RefineableElement::max_integrity_tolerance())
4756  {
4757  std::ostringstream error_stream;
4758  error_stream << "Mesh refined: Max. error in integrity check: "
4759  << max_error << " is too big"
4760  << "\ni.e. bigger than RefineableElement::"
4761  << "max_integrity_tolerance()="
4762  << RefineableElement::max_integrity_tolerance()
4763  << std::endl;
4764 
4765  std::ofstream some_file;
4766  some_file.open("ProblemMesh.dat");
4767  for (unsigned long n = 0; n < n_node; n++)
4768  {
4769  // Get the pointer to the node
4770  Node* nod_pt = this->node_pt(n);
4771  // Get the dimension
4772  unsigned n_dim = nod_pt->ndim();
4773  // Output the coordinates
4774  for (unsigned i = 0; i < n_dim; i++)
4775  {
4776  some_file << this->node_pt(n)->x(i) << " ";
4777  }
4778  some_file << std::endl;
4779  }
4780  some_file.close();
4781 
4782  error_stream << "Documented problem mesh in ProblemMesh.dat"
4783  << std::endl;
4784 
4785  throw OomphLibError(
4786  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4787  }
4788  else
4789  {
4790  oomph_info << "Mesh refined: Max. error in integrity check: "
4791  << max_error << " is OK" << std::endl;
4792  oomph_info
4793  << "i.e. less than RefineableElement::max_integrity_tolerance()="
4794  << RefineableElement::max_integrity_tolerance() << "\n"
4795  << std::endl;
4796  }
4797 
4798 
4799  if (Global_timings::Doc_comprehensive_timings)
4800  {
4801  t_end = TimingHelpers::timer();
4802  oomph_info << "Time for (paranoid only) checking of integrity: "
4803  << t_end - t_start << std::endl;
4804  t_start = TimingHelpers::timer();
4805  }
4806 
4807 #endif
4808 
4809  // Loop over all elements other than the final level and deactivate the
4810  // objects, essentially set the pointer that point to nodes that are
4811  // about to be deleted to NULL. This must take place here because nodes
4812  // addressed by elements that are dead but still living in the tree might
4813  // have been made obsolete in the last round of refinement
4814  //(Not strictly required, as tree structure has not changed, but does no
4815  // harm)
4816  for (unsigned long e = 0; e < this->forest_pt()->ntree(); e++)
4817  {
4818  this->forest_pt()->tree_pt(e)->traverse_all(&Tree::deactivate_object);
4819  }
4820 
4821  // Now we can prune the dead nodes from the mesh.
4822  Vector<Node*> deleted_node_pt = this->prune_dead_nodes();
4823 
4824  if (Global_timings::Doc_comprehensive_timings)
4825  {
4826  t_end = TimingHelpers::timer();
4827  oomph_info << "Time for deactivating objects and pruning nodes: "
4828  << t_end - t_start << std::endl;
4829  t_start = TimingHelpers::timer();
4830  }
4831 
4832  // Finally: Reorder the nodes within the mesh's node vector
4833  // to establish a standard ordering regardless of the sequence
4834  // of mesh refinements -- this is required to allow dump/restart
4835  // on refined meshes
4836  this->reorder_nodes();
4837 
4838  if (Global_timings::Doc_comprehensive_timings)
4839  {
4840  t_end = TimingHelpers::timer();
4841  oomph_info << "Time for reordering " << nnode()
4842  << " nodes: " << t_end - t_start << std::endl;
4843  t_start = TimingHelpers::timer();
4844  }
4845 
4846  // Now we can correct the nodes on boundaries that were hanging that
4847  // are no longer hanging
4848  // Only bother if we have more than two dimensions
4849  if (mesh_dim > 2)
4850  {
4851  // Loop over the boundaries
4852  for (unsigned b = 0; b < n_boundary; b++)
4853  {
4854  // Remove deleted nodes from the set
4855  unsigned n_del = deleted_node_pt.size();
4856  for (unsigned j = 0; j < n_del; j++)
4857  {
4858  hanging_nodes_on_boundary_pt[b].erase(deleted_node_pt[j]);
4859  }
4860 
4861  // If the nodes that were hanging are still hanging then remove them
4862  // from the set (note increment is not in for command for efficiencty)
4863  for (std::set<Node*>::iterator it =
4864  hanging_nodes_on_boundary_pt[b].begin();
4865  it != hanging_nodes_on_boundary_pt[b].end();)
4866  {
4867  if ((*it)->is_hanging())
4868  {
4869  hanging_nodes_on_boundary_pt[b].erase(it++);
4870  }
4871  else
4872  {
4873  ++it;
4874  }
4875  }
4876 
4877  // Are there any nodes that have changed geometric hanging status
4878  // on the boundary
4879  // The slightly painful part is that we must adjust the position
4880  // via the macro-elements which are only available through the
4881  // elements and not the nodes.
4882  if (hanging_nodes_on_boundary_pt[b].size() > 0)
4883  {
4884  // If so we loop over all elements adjacent to the boundary
4885  unsigned n_boundary_element = this->nboundary_element(b);
4886  for (unsigned e = 0; e < n_boundary_element; ++e)
4887  {
4888  // Get a pointer to the element
4889  FiniteElement* el_pt = this->boundary_element_pt(b, e);
4890 
4891  // Do we have a solid element
4892  SolidFiniteElement* solid_el_pt =
4893  dynamic_cast<SolidFiniteElement*>(el_pt);
4894 
4895  // Determine whether there is a macro element
4896  bool macro_present = (el_pt->macro_elem_pt() != 0);
4897  // Or a solid macro element
4898  if (solid_el_pt != 0)
4899  {
4900  macro_present |= (solid_el_pt->undeformed_macro_elem_pt() != 0);
4901  }
4902 
4903  // Only bother to do anything if there is a macro element
4904  // or undeformed macro element in a SolidElement
4905  if (macro_present)
4906  {
4907  // Loop over the nodes
4908  // ALH: (could optimise to only loop over
4909  // node associated with the boundary with more effort)
4910  unsigned n_el_node = el_pt->nnode();
4911  for (unsigned n = 0; n < n_el_node; n++)
4912  {
4913  // Cache pointer to the node
4914  Node* nod_pt = el_pt->node_pt(n);
4915  if (nod_pt->is_on_boundary(b))
4916  {
4917  // Is the Node in our set
4918  std::set<Node*>::iterator it =
4919  hanging_nodes_on_boundary_pt[b].find(nod_pt);
4920  // If we have found the Node then update the position
4921  // to be consistent with the macro-element representation
4922  if (it != hanging_nodes_on_boundary_pt[b].end())
4923  {
4924  // Specialise local and global coordinates to 3D
4925  // because there is only a problem in 3D.
4926  Vector<double> s(3), x(3);
4927  // Find the local coordinate of the ndoe
4928  el_pt->local_coordinate_of_node(n, s);
4929  // Find the number of time history values
4930  const unsigned ntstorage = nod_pt->ntstorage();
4931 
4932  // Do we have a solid node
4933  SolidNode* solid_node_pt =
4934  dynamic_cast<SolidNode*>(nod_pt);
4935  if (solid_node_pt)
4936  {
4937  // Assign Lagrangian coordinates from undeformed
4938  // macro element (if it has one -- get_x_and_xi()
4939  // does "the right thing" anyway. Leave actual
4940  // nodal positions alone -- we're doing a solid
4941  // mechanics problem and once we're going
4942  // the nodal positions are always computed, never
4943  // (automatically) reset to macro-element based
4944  // positions; not even on pinned boundaries
4945  // because the user may have other ideas about where
4946  // these should go -- pinning means "don't touch the
4947  // value", not "leave where the macro-element thinks
4948  // it should be"
4949  Vector<double> x_fe(3), xi(3), xi_fe(3);
4950  solid_el_pt->get_x_and_xi(s, x_fe, x, xi_fe, xi);
4951  for (unsigned i = 0; i < 3; i++)
4952  {
4953  solid_node_pt->xi(i) = xi[i];
4954  }
4955  }
4956  else
4957  {
4958  // Set position and history values from the
4959  // macro-element representation
4960  for (unsigned t = 0; t < ntstorage; t++)
4961  {
4962  // Get the history value from the macro element
4963  el_pt->get_x(t, s, x);
4964 
4965  // Set the coordinate to that of the macroelement
4966  // representation
4967  for (unsigned i = 0; i < 3; i++)
4968  {
4969  nod_pt->x(t, i) = x[i];
4970  }
4971  }
4972  } // End of non-solid node case
4973 
4974  // Now remove the node from the list
4975  hanging_nodes_on_boundary_pt[b].erase(it);
4976  // If there are no Nodes left then exit the loops
4977  if (hanging_nodes_on_boundary_pt[b].size() == 0)
4978  {
4979  e = n_boundary_element;
4980  break;
4981  }
4982  }
4983  }
4984  }
4985  } // End of macro element case
4986  }
4987  }
4988  }
4989  } // End of case when we have fixed nodal positions
4990 
4991  // Final doc
4992  //-----------
4993  if (doc_info.is_doc_enabled())
4994  {
4995  // Doc the boundary conditions ('0' for non-existent, '1' for free,
4996  //----------------------------------------------------------------
4997  // '2' for pinned -- ideal for tecplot scatter sizing.
4998  //----------------------------------------------------
4999  // num_tree_nodes=tree_nodes_pt.size();
5000 
5001  // Determine maximum number of values at any node in this type of
5002  // element
5003  RefineableElement* el_pt = tree_nodes_pt[0]->object_pt();
5004  // Initalise max_nval
5005  unsigned max_nval = 0;
5006  for (unsigned n = 0; n < el_pt->nnode(); n++)
5007  {
5008  if (el_pt->node_pt(n)->nvalue() > max_nval)
5009  {
5010  max_nval = el_pt->node_pt(n)->nvalue();
5011  }
5012  }
5013 
5014  // Open the output file
5015  std::ostringstream fullname;
5016  std::ofstream bcs_file;
5017  fullname.str("");
5018  fullname << doc_info.directory() << "/bcs" << doc_info.number()
5019  << ".dat";
5020  bcs_file.open(fullname.str().c_str());
5021 
5022  // Loop over elements
5023  for (unsigned long e = 0; e < num_tree_nodes; e++)
5024  {
5025  el_pt = tree_nodes_pt[e]->object_pt();
5026  // Loop over nodes in element
5027  unsigned n_nod = el_pt->nnode();
5028  for (unsigned n = 0; n < n_nod; n++)
5029  {
5030  // Get pointer to the node
5031  Node* nod_pt = el_pt->node_pt(n);
5032  // Find the dimension of the node
5033  unsigned n_dim = nod_pt->ndim();
5034  // Write the nodal coordinates to the file
5035  for (unsigned i = 0; i < n_dim; i++)
5036  {
5037  bcs_file << nod_pt->x(i) << " ";
5038  }
5039 
5040  // Loop over all values in this element
5041  for (unsigned i = 0; i < max_nval; i++)
5042  {
5043  // Value exists at this node:
5044  if (i < nod_pt->nvalue())
5045  {
5046  bcs_file << " " << 1 + nod_pt->is_pinned(i);
5047  }
5048  // ...if not just dump out a zero
5049  else
5050  {
5051  bcs_file << " 0 ";
5052  }
5053  }
5054  bcs_file << std::endl;
5055  }
5056  }
5057  bcs_file.close();
5058 
5059  // Doc all nodes
5060  //---------------
5061  std::ofstream all_nodes_file;
5062  fullname.str("");
5063  fullname << doc_info.directory() << "/all_nodes" << doc_info.number()
5064  << ".dat";
5065  all_nodes_file.open(fullname.str().c_str());
5066 
5067  all_nodes_file << "ZONE \n";
5068 
5069  // Need to recompute the number of nodes since it may have
5070  // changed during mesh refinement/unrefinement
5071  n_node = this->nnode();
5072  for (unsigned long n = 0; n < n_node; n++)
5073  {
5074  Node* nod_pt = this->node_pt(n);
5075  unsigned n_dim = nod_pt->ndim();
5076  for (unsigned i = 0; i < n_dim; i++)
5077  {
5078  all_nodes_file << this->node_pt(n)->x(i) << " ";
5079  }
5080  all_nodes_file << std::endl;
5081  }
5082 
5083  all_nodes_file.close();
5084 
5085 
5086  // Doc all hanging nodes:
5087  //-----------------------
5088  std::ofstream some_file;
5089  fullname.str("");
5090  fullname << doc_info.directory() << "/all_hangnodes"
5091  << doc_info.number() << ".dat";
5092  some_file.open(fullname.str().c_str());
5093  for (unsigned long n = 0; n < n_node; n++)
5094  {
5095  Node* nod_pt = this->node_pt(n);
5096 
5097  if (nod_pt->is_hanging())
5098  {
5099  unsigned n_dim = nod_pt->ndim();
5100  for (unsigned i = 0; i < n_dim; i++)
5101  {
5102  some_file << nod_pt->x(i) << " ";
5103  }
5104 
5105  // ALH: Added this to stop Solid problems seg-faulting
5106  if (this->node_pt(n)->nvalue() > 0)
5107  {
5108  some_file << " " << nod_pt->raw_value(0);
5109  }
5110  some_file << std::endl;
5111  }
5112  }
5113  some_file.close();
5114 
5115  // Doc all hanging nodes and their masters
5116  // View with QHangingNodesWithMasters.mcr
5117  fullname.str("");
5118  fullname << doc_info.directory() << "/geometric_hangnodes_withmasters"
5119  << doc_info.number() << ".dat";
5120  some_file.open(fullname.str().c_str());
5121  for (unsigned long n = 0; n < n_node; n++)
5122  {
5123  Node* nod_pt = this->node_pt(n);
5124  if (nod_pt->is_hanging())
5125  {
5126  unsigned n_dim = nod_pt->ndim();
5127  unsigned nmaster = nod_pt->hanging_pt()->nmaster();
5128  some_file << "ZONE I=" << nmaster + 1 << std::endl;
5129  for (unsigned i = 0; i < n_dim; i++)
5130  {
5131  some_file << nod_pt->x(i) << " ";
5132  }
5133  some_file << " 2 " << std::endl;
5134 
5135  for (unsigned imaster = 0; imaster < nmaster; imaster++)
5136  {
5137  Node* master_nod_pt =
5138  nod_pt->hanging_pt()->master_node_pt(imaster);
5139  unsigned n_dim = master_nod_pt->ndim();
5140  for (unsigned i = 0; i < n_dim; i++)
5141  {
5142  some_file << master_nod_pt->x(i) << " ";
5143  }
5144  some_file << " 1 " << std::endl;
5145  }
5146  }
5147  }
5148  some_file.close();
5149 
5150  // Doc all hanging nodes and their masters
5151  // View with QHangingNodesWithMasters.mcr
5152  for (unsigned i = 0; i < ncont_interpolated_values; i++)
5153  {
5154  fullname.str("");
5155  fullname << doc_info.directory()
5156  << "/nonstandard_hangnodes_withmasters" << i << "_"
5157  << doc_info.number() << ".dat";
5158  some_file.open(fullname.str().c_str());
5159  unsigned n_nod = this->nnode();
5160  for (unsigned long n = 0; n < n_nod; n++)
5161  {
5162  Node* nod_pt = this->node_pt(n);
5163  if (nod_pt->is_hanging(i))
5164  {
5165  if (nod_pt->hanging_pt(i) != nod_pt->hanging_pt())
5166  {
5167  unsigned nmaster = nod_pt->hanging_pt(i)->nmaster();
5168  some_file << "ZONE I=" << nmaster + 1 << std::endl;
5169  unsigned n_dim = nod_pt->ndim();
5170  for (unsigned j = 0; j < n_dim; j++)
5171  {
5172  some_file << nod_pt->x(j) << " ";
5173  }
5174  some_file << " 2 " << std::endl;
5175  for (unsigned imaster = 0; imaster < nmaster; imaster++)
5176  {
5177  Node* master_nod_pt =
5178  nod_pt->hanging_pt(i)->master_node_pt(imaster);
5179  unsigned n_dim = master_nod_pt->ndim();
5180  for (unsigned j = 0; j < n_dim; j++)
5181  {
5182  // some_file << master_nod_pt->x(i) << " ";
5183  }
5184  some_file << " 1 " << std::endl;
5185  }
5186  }
5187  }
5188  }
5189  some_file.close();
5190  }
5191 
5192  } // End of documentation
5193  } // End if (this->nelement()>0)
5194 
5196  // std::cout << "p_adapt_mesh(): Checking stuff works!" << std::endl;
5197  // for(unsigned j=0; j<this->nnode(); j++)
5198  // {
5199  // MacroElementNodeUpdateNode* macro_nod_pt =
5200  // dynamic_cast<MacroElementNodeUpdateNode*>(this->node_pt(j));
5201  // if(macro_nod_pt!=0)
5202  // {
5203  // bool big_problem = true;
5204  // std::cout << "Node " << macro_nod_pt << " at [ " << macro_nod_pt->x(0)
5205  // << ", " << macro_nod_pt->x(1) << " ]" << std::endl; FiniteElement*
5206  // up_el_pt =
5207  // dynamic_cast<FiniteElement*>(macro_nod_pt->node_update_element_pt());
5208  // for(unsigned l=0; l<up_el_pt->nnode(); l++)
5209  // {
5210  // if(up_el_pt->node_pt(l)==macro_nod_pt)
5211  // {
5212  // big_problem = false;
5213  // break;
5214  // }
5215  // }
5216  // if(big_problem)
5217  // {
5218  // std::cout << " This node doesn't exist in it's update element!" <<
5219  // std::endl;
5220  // }
5221  // }
5222  // }
5223 
5224 #ifdef OOMPH_HAS_MPI
5225 
5226  // Now (re-)classify halo and haloed nodes and synchronise hanging
5227  // nodes
5228  if (this->is_mesh_distributed())
5229  {
5230  classify_halo_and_haloed_nodes(doc_info, doc_info.is_doc_enabled());
5231  }
5232 
5233 #endif
5234  }
oomph::TreeBasedRefineableMeshBase::p_refine_elements_if_required
virtual void p_refine_elements_if_required()=0

References b, oomph::Mesh::boundary_element_pt(), oomph::TreeForest::check_all_neighbours(), oomph::TreeForest::close_hanging_node_files(), complete_hanging_nodes(), oomph::Tree::deactivate_object(), oomph::Mesh::delete_all_external_storage(), oomph::PRefineableElement::deselect_for_p_refinement(), oomph::PRefineableElement::deselect_for_p_unrefinement(), oomph::FiniteElement::dim(), oomph::DocInfo::directory(), oomph::Global_timings::Doc_comprehensive_timings, oomph::RefineableMeshBase::doc_info(), e(), oomph::Mesh::element_pt(), Eigen::placeholders::end, oomph::Mesh::finite_element_pt(), forest_pt(), oomph::MacroElementNodeUpdateElementBase::geom_object_pt(), oomph::Node::get_boundaries_pt(), oomph::FiniteElement::get_x(), oomph::SolidFiniteElement::get_x_and_xi(), oomph::Node::hanging_pt(), i, oomph::DocInfo::is_doc_enabled(), oomph::Node::is_hanging(), oomph::Mesh::is_mesh_distributed(), oomph::Node::is_on_boundary(), oomph::Data::is_pinned(), j, oomph::SolidNode::lagrangian_position(), oomph::FiniteElement::local_coordinate_of_node(), oomph::Mesh::Lookup_for_elements_next_boundary_is_setup, oomph::FiniteElement::macro_elem_pt(), oomph::HangInfo::master_node_pt(), oomph::RefineableMeshBase::max_error(), oomph::RefineableElement::max_integrity_tolerance(), n, oomph::Mesh::nboundary(), oomph::Mesh::nboundary_element(), oomph::Node::ndim(), oomph::Mesh::nelement(), oomph::SolidNode::nlagrangian(), oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::Mesh::nnode(), oomph::FiniteElement::node_pt(), oomph::Mesh::node_pt(), oomph::AlgebraicNode::node_update(), oomph::TreeForest::ntree(), oomph::Data::ntstorage(), oomph::DocInfo::number(), oomph::Data::nvalue(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::TreeForest::open_hanging_node_files(), p_refine_elements_if_required(), oomph::Node::position(), oomph::Mesh::prune_dead_nodes(), oomph::Node::raw_value(), oomph::Mesh::reorder_nodes(), s, oomph::MacroElementNodeUpdateElementBase::set_node_update_info(), oomph::Node::set_non_obsolete(), oomph::Node::set_nonhanging(), oomph::Node::set_obsolete(), oomph::Data::set_value(), oomph::Mesh::setup_boundary_element_info(), size, oomph::TreeForest::stick_all_tree_nodes_into_vector(), oomph::TreeForest::stick_leaves_into_vector(), plotPSD::t, oomph::TimingHelpers::timer(), oomph::Tree::traverse_all(), oomph::TreeForest::tree_pt(), oomph::SolidFiniteElement::undeformed_macro_elem_pt(), oomph::Node::value(), plotDoE::x, oomph::Node::x(), and oomph::SolidNode::xi().

◆ p_refine_elements_if_required()

virtual void oomph::TreeBasedRefineableMeshBase::p_refine_elements_if_required ( )
protectedpure virtual

p-refine all the elements in the mesh if required. This template free interface will be overloaded in RefineableMesh<ELEMENT> so that any temporary copies of the element that are created will be of the correct type.

Implemented in oomph::TreeBasedRefineableMesh< ELEMENT >.

Referenced by p_adapt_mesh().

◆ p_refine_selected_elements() [1/2]

void oomph::TreeBasedRefineableMeshBase::p_refine_selected_elements ( const Vector< PRefineableElement * > &  elements_to_be_refined_pt)

p-refine mesh by refining the elements identified by their pointers.

p-refine mesh by refining the elements identified by their pointers.

5306  {
5307 #ifdef OOMPH_HAS_MPI
5308  if (this->is_mesh_distributed())
5309  {
5310  std::ostringstream warn_stream;
5311  warn_stream << "You are attempting to refine selected elements of a "
5312  << std::endl
5313  << "distributed mesh. This may have undesired effects."
5314  << std::endl;
5315 
5316  OomphLibWarning(warn_stream.str(),
5317  "TreeBasedRefineableMeshBase::refine_selected_elements()",
5318  OOMPH_EXCEPTION_LOCATION);
5319  }
5320 #endif
5321 
5322  // Select elements for refinement
5323  unsigned long nref = elements_to_be_refined_pt.size();
5324  for (unsigned long e = 0; e < nref; e++)
5325  {
5326  elements_to_be_refined_pt[e]->select_for_p_refinement();
5327  }
5328 
5329  // Do the actual mesh adaptation
5330  p_adapt_mesh();
5331  }

References e(), oomph::Mesh::is_mesh_distributed(), OOMPH_EXCEPTION_LOCATION, and p_adapt_mesh().

◆ p_refine_selected_elements() [2/2]

void oomph::TreeBasedRefineableMeshBase::p_refine_selected_elements ( const Vector< unsigned > &  elements_to_be_refined)

p-refine mesh by refining the elements identified by their numbers.

p-refine mesh by refining the elements identified by their numbers.

5266  {
5267 #ifdef OOMPH_HAS_MPI
5268  if (this->is_mesh_distributed())
5269  {
5270  std::ostringstream warn_stream;
5271  warn_stream << "You are attempting to refine selected elements of a "
5272  << std::endl
5273  << "distributed mesh. This may have undesired effects."
5274  << std::endl;
5275 
5276  OomphLibWarning(warn_stream.str(),
5277  "TreeBasedRefineableMeshBase::refine_selected_elements()",
5278  OOMPH_EXCEPTION_LOCATION);
5279  }
5280 #endif
5281 
5282  // Select elements for refinement
5283  unsigned long nref = elements_to_be_refined.size();
5284  for (unsigned long e = 0; e < nref; e++)
5285  {
5286  // Get pointer to p-refineable element
5287  PRefineableElement* el_pt = dynamic_cast<PRefineableElement*>(
5288  this->element_pt(elements_to_be_refined[e]));
5289  // Mark for p-refinement if possible. If not then p_adapt_mesh() will
5290  // report the error.
5291  if (el_pt != 0)
5292  {
5293  el_pt->select_for_p_refinement();
5294  }
5295  }
5296 
5297  // Do the actual mesh adaptation
5298  p_adapt_mesh();
5299  }

References e(), oomph::Mesh::element_pt(), oomph::Mesh::is_mesh_distributed(), OOMPH_EXCEPTION_LOCATION, p_adapt_mesh(), and oomph::PRefineableElement::select_for_p_refinement().

◆ p_refine_uniformly() [1/2]

void oomph::TreeBasedRefineableMeshBase::p_refine_uniformly ( )
inlinevirtual

p-refine mesh uniformly

Reimplemented from oomph::RefineableMeshBase.

444  {
445  RefineableMeshBase::p_refine_uniformly();
446  }
oomph::RefineableMeshBase::p_refine_uniformly
virtual void p_refine_uniformly()
p-refine mesh uniformly
Definition: refineable_mesh.h:294

References oomph::RefineableMeshBase::p_refine_uniformly().

◆ p_refine_uniformly() [2/2]

void oomph::TreeBasedRefineableMeshBase::p_refine_uniformly ( DocInfo &  doc_info)
virtual

p-refine mesh uniformly and doc process

p-refine mesh uniformly

Reimplemented from oomph::RefineableMeshBase.

1791  {
1792  // Select all elements for refinement
1793  unsigned long Nelement = this->nelement();
1794  for (unsigned long e = 0; e < Nelement; e++)
1795  {
1796  // Get pointer to p-refineable element
1797  PRefineableElement* el_pt =
1798  dynamic_cast<PRefineableElement*>(this->element_pt(e));
1799  // Mark for p-refinement if possible. If not then p_adapt_mesh() will
1800  // report the error.
1801  if (el_pt != 0)
1802  {
1803  el_pt->select_for_p_refinement();
1804  }
1805  }
1806 
1807  // Do the actual mesh adaptation
1808  p_adapt_mesh(doc_info);
1809  }

References oomph::RefineableMeshBase::doc_info(), e(), oomph::Mesh::element_pt(), oomph::Mesh::nelement(), p_adapt_mesh(), and oomph::PRefineableElement::select_for_p_refinement().

◆ p_unrefine_uniformly()

void oomph::TreeBasedRefineableMeshBase::p_unrefine_uniformly ( DocInfo &  doc_info)

p-unrefine mesh uniformly

p-unrefine mesh uniformly Unlike in h-refinement, we can simply p-unrefine each element in the mesh

5241  {
5242  // Select all elements for unrefinement
5243  unsigned long Nelement = this->nelement();
5244  for (unsigned long e = 0; e < Nelement; e++)
5245  {
5246  // Get pointer to p-refineable element
5247  PRefineableElement* el_pt =
5248  dynamic_cast<PRefineableElement*>(this->element_pt(e));
5249  // Mark for p-refinement if possible. If not then p_adapt_mesh() will
5250  // report the error.
5251  if (el_pt != 0)
5252  {
5253  el_pt->select_for_p_unrefinement();
5254  }
5255  }
5256 
5257  // Do the actual mesh adaptation
5258  p_adapt_mesh(doc_info);
5259  }

References oomph::RefineableMeshBase::doc_info(), e(), oomph::Mesh::element_pt(), oomph::Mesh::nelement(), p_adapt_mesh(), and oomph::PRefineableElement::select_for_p_unrefinement().

◆ read_refinement()

void oomph::TreeBasedRefineableMeshBase::read_refinement ( std::ifstream &  restart_file,
Vector< Vector< unsigned >> &  to_be_refined 
)
virtual

Read refinement pattern to allow for rebuild.

Read refinement pattern to allow for rebuild

253  {
254  std::string input_string;
255 
256  // Read max refinement level:
257 
258  // Read line up to termination sign
259  getline(restart_file, input_string, '#');
260 
261  // Ignore rest of line
262  restart_file.ignore(80, '\n');
263 
264  // Convert
265  unsigned max_level = std::atoi(input_string.c_str());
266 
267  // Assign storage for refinement pattern
268  to_be_refined.resize(max_level);
269 
270  // Read the number of the elements that need to be refined at different
271  // levels
272  for (unsigned l = 0; l < max_level; l++)
273  {
274  // Read line up to termination sign
275  getline(restart_file, input_string, '#');
276 
277  // Ignore rest of line
278  restart_file.ignore(80, '\n');
279 
280  // Convert
281  unsigned n_to_be_refined = atoi(input_string.c_str());
282  ;
283 
284  // Assign storage
285  to_be_refined[l].resize(n_to_be_refined);
286 
287  // Read numbers of the elements that need to be refined
288  for (unsigned i = 0; i < n_to_be_refined; i++)
289  {
290  restart_file >> to_be_refined[l][i];
291  }
292  }
293  }
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286

References i, and oomph::Global_string_for_annotation::string().

Referenced by refine().

◆ refine()

void oomph::TreeBasedRefineableMeshBase::refine ( std::ifstream &  restart_file)
virtual

Refine mesh according to refinement pattern in restart file.

Refine base mesh according to refinement pattern in restart file.

202  {
203  // Assign storage for refinement pattern
204  Vector<Vector<unsigned>> to_be_refined;
205 
206  // Read refinement pattern
207  read_refinement(restart_file, to_be_refined);
208 
209  // Refine
210  refine_base_mesh(to_be_refined);
211  }
oomph::TreeBasedRefineableMeshBase::refine_base_mesh
void refine_base_mesh(Vector< Vector< unsigned >> &to_be_refined)
Refine base mesh according to specified refinement pattern.
Definition: refineable_mesh.cc:137
oomph::TreeBasedRefineableMeshBase::read_refinement
virtual void read_refinement(std::ifstream &restart_file, Vector< Vector< unsigned >> &to_be_refined)
Read refinement pattern to allow for rebuild.
Definition: refineable_mesh.cc:251

References read_refinement(), and refine_base_mesh().

◆ refine_as_in_reference_mesh()

void oomph::TreeBasedRefineableMeshBase::refine_as_in_reference_mesh ( TreeBasedRefineableMeshBase *const &  ref_mesh_pt)
virtual

Refine mesh once so that its topology etc becomes that of the (finer!) reference mesh – if possible! Useful for meshes in multigrid hierarchies. If the meshes are too different and the conversion cannot be performed, the code dies (provided PARANOID is enabled).

1949  {
1950  oomph_info << "WARNING : This has not been checked comprehensively yet"
1951  << std::endl
1952  << "Check it and remove this break " << std::endl;
1953  pause("Yes really pause");
1954 
1955 #ifdef PARANOID
1956  // The max. refinement levels of the two meshes need to differ
1957  // by one, otherwise what we're doing here doesn't make sense.
1958  unsigned my_min, my_max;
1959  get_refinement_levels(my_min, my_max);
1960 
1961  unsigned ref_min, ref_max;
1962  ref_mesh_pt->get_refinement_levels(ref_min, ref_max);
1963 
1964  if (ref_max != my_max + 1)
1965  {
1966  std::ostringstream error_stream;
1967  error_stream
1968  << "Meshes definitely don't differ by one refinement level \n"
1969  << "max. refinement levels: " << ref_max << " " << my_max << std::endl;
1970 
1971  throw OomphLibError(
1972  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1973  }
1974 #endif
1975 
1976  // Vector storing the elements of the uniformly unrefined mesh
1977  Vector<Tree*> coarse_elements_pt;
1978 
1979  // Map storing which father elements have already been added to coarse mesh
1980  // (Default return is 0).
1981  std::map<Tree*, unsigned> father_element_included;
1982 
1983  // Extract active elements (=leaf nodes in the quadtree) from reference
1984  // mesh.
1985  Vector<Tree*> leaf_nodes_pt;
1986  ref_mesh_pt->forest_pt()->stick_leaves_into_vector(leaf_nodes_pt);
1987 
1988  // Loop over all elements (in their quadtree impersonation) and
1989  // check if their fathers's sons are all leaves too:
1990  unsigned nelem = leaf_nodes_pt.size();
1991  for (unsigned e = 0; e < nelem; e++)
1992  {
1993  // Pointer to leaf node
1994  Tree* leaf_pt = leaf_nodes_pt[e];
1995 
1996  // Get pointer to father:
1997  Tree* father_pt = leaf_pt->father_pt();
1998 
1999  // If we don't have a father we're at the root level in which
2000  // case this element can't be unrefined.
2001  if (0 == father_pt)
2002  {
2003  coarse_elements_pt.push_back(leaf_pt);
2004  }
2005  else
2006  {
2007  // Loop over the father's sons to check if they're
2008  // all non-leafs, i.e. if they can be unrefined
2009  bool can_unrefine = true;
2010  unsigned n_sons = father_pt->nsons();
2011  for (unsigned i = 0; i < n_sons; i++)
2012  {
2013  // If (at least) one of the sons is not a leaf, we can't unrefine
2014  if (!father_pt->son_pt(i)->is_leaf()) can_unrefine = false;
2015  }
2016 
2017  // If we can unrefine, the father element will be
2018  // an element in the coarse mesh, the sons won't
2019  if (can_unrefine)
2020  {
2021  if (father_element_included[father_pt] == 0)
2022  {
2023  coarse_elements_pt.push_back(father_pt);
2024  father_element_included[father_pt] = 1;
2025  }
2026  }
2027  // Son will still be there on the coarse mesh
2028  else
2029  {
2030  coarse_elements_pt.push_back(leaf_pt);
2031  }
2032  }
2033  }
2034 
2035  // Number of elements in ref mesh if it was unrefined uniformly:
2036  unsigned nel_coarse = coarse_elements_pt.size();
2037 
2038 
2039 #ifdef PARANOID
2040  bool stop_it = false;
2041  // The numbers had better match otherwise we might as well stop now...
2042  if (nel_coarse != this->nelement())
2043  {
2044  oomph_info << "Number of elements in uniformly unrefined reference mesh: "
2045  << nel_coarse << std::endl;
2046  oomph_info << "Number of elements in 'this' mesh: " << nel_coarse
2047  << std::endl;
2048  oomph_info << "don't match" << std::endl;
2049  stop_it = true;
2050  }
2051 #endif
2052 
2053  // Now loop over all elements in uniformly coarsened reference mesh
2054  // and check if add the number of any element that was created
2055  // by having had its sons merged to the vector of elements that
2056  // need to get refined if we go the other way
2057  Vector<unsigned> elements_to_be_refined;
2058  for (unsigned i = 0; i < nel_coarse; i++)
2059  {
2060  if (father_element_included[coarse_elements_pt[i]] == 1)
2061  {
2062  elements_to_be_refined.push_back(i);
2063  }
2064  }
2065 
2066 
2067 #ifdef PARANOID
2068  // Doc troublesome meshes:
2069  if (stop_it)
2070  {
2071  std::ofstream some_file;
2072  some_file.open("orig_mesh.dat");
2073  this->output(some_file);
2074  some_file.close();
2075  oomph_info << "Documented original ('this')mesh in orig_mesh.dat"
2076  << std::endl;
2077  }
2078 #endif
2079 
2080 
2081  // Now refine precisely these elements in "this" mesh.
2082  refine_selected_elements(elements_to_be_refined);
2083 
2084 
2085 #ifdef PARANOID
2086 
2087  // Check if the nodal positions of all element's nodes agree
2088  // in the two fine meshes:
2089  double tol = 1.0e-5;
2090  for (unsigned e = 0; e < nelem; e++)
2091  {
2092  // Get elements
2093  FiniteElement* ref_el_pt = ref_mesh_pt->finite_element_pt(e);
2094  FiniteElement* el_pt = this->finite_element_pt(e);
2095 
2096  // Loop over nodes
2097  unsigned nnod = ref_el_pt->nnode();
2098  for (unsigned j = 0; j < nnod; j++)
2099  {
2100  // Get nodes
2101  Node* ref_node_pt = ref_el_pt->node_pt(j);
2102  Node* node_pt = el_pt->node_pt(j);
2103 
2104  // Check error in position
2105  double error = 0.0;
2106  unsigned ndim = node_pt->ndim();
2107  for (unsigned i = 0; i < ndim; i++)
2108  {
2109  error += pow(node_pt->x(i) - ref_node_pt->x(i), 2);
2110  }
2111  error = sqrt(error);
2112 
2113  if (error > tol)
2114  {
2115  oomph_info << "Error in nodal position of node " << j << ": " << error
2116  << " [tol=" << tol << "]" << std::endl;
2117  stop_it = true;
2118  }
2119  }
2120  }
2121 
2122  // Do we have a death wish?
2123  if (stop_it)
2124  {
2125  // Doc troublesome meshes:
2126  std::ofstream some_file;
2127  some_file.open("refined_mesh.dat");
2128  this->output(some_file);
2129  some_file.close();
2130 
2131  some_file.open("finer_mesh.dat");
2132  ref_mesh_pt->output(some_file);
2133  some_file.close();
2134 
2135  throw OomphLibError(
2136  "Bailing out. Doced refined_mesh.dat finer_mesh.dat\n",
2137  OOMPH_CURRENT_FUNCTION,
2138  OOMPH_EXCEPTION_LOCATION);
2139  }
2140 
2141 #endif
2142  }
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
oomph::Mesh::output
void output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:2027
oomph::Node::ndim
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
oomph::TreeBasedRefineableMeshBase::refine_selected_elements
virtual void refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
Definition: refineable_mesh.cc:1816
oomph::TreeBasedRefineableMeshBase::get_refinement_levels
virtual void get_refinement_levels(unsigned &min_refinement_level, unsigned &max_refinement_level)
Get max/min refinement levels in mesh.
Definition: refineable_mesh.cc:830
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
calibrate.error
int error
Definition: calibrate.py:297
oomph::pause
void pause(std::string message)
Pause and display message.
Definition: oomph_utilities.cc:1265

References e(), calibrate::error, oomph::Tree::father_pt(), oomph::Mesh::finite_element_pt(), forest_pt(), get_refinement_levels(), i, oomph::Tree::is_leaf(), j, oomph::Node::ndim(), oomph::Mesh::nelement(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), oomph::Mesh::node_pt(), oomph::Tree::nsons(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::Mesh::output(), oomph::pause(), Eigen::bfloat16_impl::pow(), refine_selected_elements(), oomph::Tree::son_pt(), sqrt(), oomph::TreeForest::stick_leaves_into_vector(), and oomph::Node::x().

◆ refine_base_mesh()

void oomph::TreeBasedRefineableMeshBase::refine_base_mesh ( Vector< Vector< unsigned >> &  to_be_refined)

Refine base mesh according to specified refinement pattern.

Refine original, unrefined mesh according to specified refinement pattern (relative to original, unrefined mesh).

139  {
140  // Get mesh back to unrefined state
141  unsigned my_max, my_min;
142  get_refinement_levels(my_min, my_max);
143 
144  // Max refinement level:
145  unsigned my_max_level = to_be_refined.size();
146 
147  unsigned global_max = 0;
148  unsigned global_max_level = 0;
149  Vector<unsigned> data(2, 0);
150  data[0] = my_max;
151  data[1] = my_max_level;
152  Vector<unsigned> global_data(2, 0);
153 #ifdef OOMPH_HAS_MPI
154  if (this->is_mesh_distributed())
155  {
156  MPI_Allreduce(&data[0],
157  &global_data[0],
158  2,
159  MPI_UNSIGNED,
160  MPI_MAX,
161  Comm_pt->mpi_comm());
162  global_max = global_data[0];
163  global_max_level = global_data[1];
164  }
165  else
166 #endif
167  {
168  global_max = my_max;
169  global_max_level = my_max_level;
170  }
171 
172 
173  for (unsigned i = 0; i < global_max; i++)
174  {
175  unrefine_uniformly();
176  }
177 
178  // Do refinement steps in current mesh
179  for (unsigned l = 0; l < global_max_level; l++)
180  {
181  // Loop over elements that need to be refined at this level
182  unsigned n_to_be_refined = 0;
183  if (l < my_max_level) n_to_be_refined = to_be_refined[l].size();
184 
185  // Select relevant elements to be refined
186  for (unsigned i = 0; i < n_to_be_refined; i++)
187  {
188  dynamic_cast<RefineableElement*>(this->element_pt(to_be_refined[l][i]))
189  ->select_for_refinement();
190  }
191 
192  // Now do the actual mesh refinement
193  adapt_mesh();
194  }
195  }
data
int data[]
Definition: Map_placement_new.cpp:1
oomph::TreeBasedRefineableMeshBase::unrefine_uniformly
unsigned unrefine_uniformly()
Definition: refineable_mesh.cc:2150

References adapt_mesh(), data, oomph::Mesh::element_pt(), get_refinement_levels(), i, oomph::Mesh::is_mesh_distributed(), and unrefine_uniformly().

Referenced by PolarNSProblem< ELEMENT >::output_streamfunction(), refine(), refine_base_mesh_as_in_reference_mesh(), and refine_base_mesh_as_in_reference_mesh_minus_one().

◆ refine_base_mesh_as_in_reference_mesh()

void oomph::TreeBasedRefineableMeshBase::refine_base_mesh_as_in_reference_mesh ( TreeBasedRefineableMeshBase *const &  ref_mesh_pt)
virtual

Refine to same degree as the reference mesh.

Refine base mesh to same degree as reference mesh (relative to original unrefined mesh).

1887  {
1888  // Assign storage for refinement pattern
1889  Vector<Vector<unsigned>> to_be_refined;
1890 
1891  // Get refinement pattern of reference mesh:
1892  ref_mesh_pt->get_refinement_pattern(to_be_refined);
1893 
1894  // Refine mesh according to given refinement pattern
1895  refine_base_mesh(to_be_refined);
1896  }

References get_refinement_pattern(), and refine_base_mesh().

◆ refine_base_mesh_as_in_reference_mesh_minus_one()

bool oomph::TreeBasedRefineableMeshBase::refine_base_mesh_as_in_reference_mesh_minus_one ( TreeBasedRefineableMeshBase *const &  ref_mesh_pt)
virtual

Refine base mesh to same degree as reference mesh minus one level of refinement (relative to original unrefined mesh). Useful function for multigrid solvers; allows the easy copy of a mesh to the level of refinement just below the current one

Refine to same degree as the reference mesh minus one. Useful function for multigrid solvers; allows the easy copy of a mesh to the level of refinement just below the current one. Returns a boolean variable which indicates if the reference mesh has not been refined at all

1909  {
1910  // Assign storage for refinement pattern
1911  Vector<Vector<unsigned>> to_be_refined;
1912 
1913  // Get refinement pattern of reference mesh:
1914  ref_mesh_pt->get_refinement_pattern(to_be_refined);
1915 
1916  // Find the length of the vector
1917  unsigned nrefinement_levels = to_be_refined.size();
1918 
1919  // If the reference mesh has not been refined a single time then
1920  // we cannot create an unrefined copy so stop here
1921  if (nrefinement_levels == 0)
1922  {
1923  return false;
1924  }
1925  // If the reference mesh has been refined at least once
1926  else
1927  {
1928  // Remove the last entry of the vector to make sure we refine to
1929  // the same level minus one
1930  to_be_refined.resize(nrefinement_levels - 1);
1931 
1932  // Refine mesh according to given refinement pattern
1933  refine_base_mesh(to_be_refined);
1934 
1935  // Indicate that it was possible to create an unrefined copy
1936  return true;
1937  }
1938  }

References get_refinement_pattern(), and refine_base_mesh().

Referenced by oomph::MGSolver< DIM >::setup_mg_hierarchy(), and oomph::HelmholtzMGPreconditioner< DIM >::setup_mg_hierarchy().

◆ refine_selected_elements() [1/2]

void oomph::TreeBasedRefineableMeshBase::refine_selected_elements ( const Vector< RefineableElement * > &  elements_to_be_refined_pt)
virtual

Refine mesh by splitting the elements identified by their pointers.

Refine mesh by splitting the elements identified by their pointers

1854  {
1855 #ifdef OOMPH_HAS_MPI
1856  if (this->is_mesh_distributed())
1857  {
1858  std::ostringstream warn_stream;
1859  warn_stream << "You are attempting to refine selected elements of a "
1860  << std::endl
1861  << "distributed mesh. This may have undesired effects."
1862  << std::endl;
1863 
1864  OomphLibWarning(warn_stream.str(),
1865  "TreeBasedRefineableMeshBase::refine_selected_elements()",
1866  OOMPH_EXCEPTION_LOCATION);
1867  }
1868 #endif
1869 
1870  // Select elements for refinement
1871  unsigned long nref = elements_to_be_refined_pt.size();
1872  for (unsigned long e = 0; e < nref; e++)
1873  {
1874  elements_to_be_refined_pt[e]->select_for_refinement();
1875  }
1876 
1877  // Do the actual mesh adaptation
1878  adapt_mesh();
1879  }

References adapt_mesh(), e(), oomph::Mesh::is_mesh_distributed(), and OOMPH_EXCEPTION_LOCATION.

◆ refine_selected_elements() [2/2]

void oomph::TreeBasedRefineableMeshBase::refine_selected_elements ( const Vector< unsigned > &  elements_to_be_refined)
virtual

Refine mesh by splitting the elements identified by their numbers.

1818  {
1819 #ifdef OOMPH_HAS_MPI
1820  if (this->is_mesh_distributed())
1821  {
1822  std::ostringstream warn_stream;
1823  warn_stream << "You are attempting to refine selected elements of a "
1824  << std::endl
1825  << "distributed mesh. This may have undesired effects."
1826  << std::endl;
1827 
1828  OomphLibWarning(warn_stream.str(),
1829  "TreeBasedRefineableMeshBase::refine_selected_elements()",
1830  OOMPH_EXCEPTION_LOCATION);
1831  }
1832 #endif
1833 
1834  // Select elements for refinement
1835  unsigned long nref = elements_to_be_refined.size();
1836  for (unsigned long e = 0; e < nref; e++)
1837  {
1838  dynamic_cast<RefineableElement*>(
1839  this->element_pt(elements_to_be_refined[e]))
1840  ->select_for_refinement();
1841  }
1842 
1843  // Do the actual mesh adaptation
1844  adapt_mesh();
1845  }

References adapt_mesh(), e(), oomph::Mesh::element_pt(), oomph::Mesh::is_mesh_distributed(), and OOMPH_EXCEPTION_LOCATION.

Referenced by FSIRingProblem::FSIRingProblem(), and refine_as_in_reference_mesh().

◆ refine_uniformly() [1/2]

void oomph::TreeBasedRefineableMeshBase::refine_uniformly ( )
inlinevirtual

Refine mesh uniformly.

Reimplemented from oomph::RefineableMeshBase.

Reimplemented in oomph::RefineableQuadFromTriangleMesh< ELEMENT >.

435  {
436  RefineableMeshBase::refine_uniformly();
437  }
oomph::RefineableMeshBase::refine_uniformly
virtual void refine_uniformly()
Refine mesh uniformly.
Definition: refineable_mesh.h:271

References oomph::RefineableMeshBase::refine_uniformly().

◆ refine_uniformly() [2/2]

void oomph::TreeBasedRefineableMeshBase::refine_uniformly ( DocInfo &  doc_info)
virtual

Refine mesh uniformly and doc process.

Refine mesh uniformly.

Implements oomph::RefineableMeshBase.

Reimplemented in oomph::RefineableQuadFromTriangleMesh< ELEMENT >.

1773  {
1774  // Select all elements for refinement
1775  unsigned long Nelement = this->nelement();
1776  for (unsigned long e = 0; e < Nelement; e++)
1777  {
1778  dynamic_cast<RefineableElement*>(this->element_pt(e))
1779  ->select_for_refinement();
1780  }
1781 
1782  // Do the actual mesh adaptation
1783  adapt_mesh(doc_info);
1784  }

References adapt_mesh(), oomph::RefineableMeshBase::doc_info(), e(), oomph::Mesh::element_pt(), and oomph::Mesh::nelement().

Referenced by ExtrudedMovingCylinderProblem< TWO_D_ELEMENT, THREE_D_ELEMENT >::create_extruded_mesh(), run_navier_stokes_outflow(), and TurekProblem< FLUID_ELEMENT, SOLID_ELEMENT >::TurekProblem().

◆ setup_tree_forest()

virtual void oomph::TreeBasedRefineableMeshBase::setup_tree_forest ( )
pure virtual

Set up the tree forest associated with the Mesh (if any)

Implemented in oomph::RefineableQuadMesh< ELEMENT >, oomph::RefineableQuadMesh< NST_ELEMENT >, oomph::RefineableQuadMesh< AD_ELEMENT >, oomph::RefineableLineMesh< ELEMENT >, and oomph::RefineableBrickMesh< ELEMENT >.

◆ split_elements_if_required()

virtual void oomph::TreeBasedRefineableMeshBase::split_elements_if_required ( )
protectedpure virtual

Split all the elements in the mesh if required. This template free interface will be overloaded in RefineableMesh<ELEMENT> so that any new elements that are created will be of the correct type.

Implemented in oomph::TreeBasedRefineableMesh< ELEMENT >.

Referenced by adapt_mesh().

◆ uniform_refinement_level_when_pruned() [1/2]

unsigned& oomph::TreeBasedRefineableMeshBase::uniform_refinement_level_when_pruned ( )
inline

Level to which the mesh was uniformly refined when it was pruned.

626  {
627  return Uniform_refinement_level_when_pruned;
628  }

References Uniform_refinement_level_when_pruned.

◆ uniform_refinement_level_when_pruned() [2/2]

unsigned oomph::TreeBasedRefineableMeshBase::uniform_refinement_level_when_pruned ( ) const
inline

Level to which the mesh was uniformly refined when it was pruned (const version)

619  {
620  return Uniform_refinement_level_when_pruned;
621  }

References Uniform_refinement_level_when_pruned.

◆ unrefine_uniformly()

unsigned oomph::TreeBasedRefineableMeshBase::unrefine_uniformly ( )
virtual

Unrefine mesh uniformly: Return 0 for success, 1 for failure (if unrefinement has reached the coarsest permitted level)

Unrefine mesh uniformly. Return 0 for success, 1 for failure (if unrefinement has reached the coarsest permitted level)

Implements oomph::RefineableMeshBase.

2151  {
2152  // We can't just select all elements for unrefinement
2153  // because they need to merge with their brothers.
2154  // --> Rather than repeating the convoluted logic of
2155  // RefineableQuadMesh<ELEMENT>::adapt(Vector<double>& elemental_error)
2156  // here (code duplication!) hack it by filling the error
2157  // vector with values that ensure unrefinement for all
2158  // elements where this is possible
2159 
2160  // Create dummy vector for elemental errors
2161  unsigned long Nelement = this->nelement();
2162  Vector<double> elemental_error(Nelement);
2163 
2164  // In order to force unrefinement, set the min permitted error to
2165  // be the default and then set the actual error to be below this.
2166  // This avoids problems when the actual min error is zero (or small)
2167  // For sanity's sake, also set the max permitted error back to the default
2168  // so that we have a max error bigger than a min error
2169  const double current_min_error = this->min_permitted_error();
2170  const double current_max_error = this->max_permitted_error();
2171 
2172  this->min_permitted_error() = 1.0e-5;
2173  this->max_permitted_error() = 1.0e-3;
2174 
2175  double error = min_permitted_error() / 100.0;
2176  for (unsigned long e = 0; e < Nelement; e++)
2177  {
2178  elemental_error[e] = error;
2179  }
2180 
2181  // Temporarily lift any restrictions on the minimum number of
2182  // elements that need to be unrefined to make it worthwhile
2183  unsigned backup = max_keep_unrefined();
2184  max_keep_unrefined() = 0;
2185 
2186  // Do the actual mesh adaptation with fake error vector
2187  adapt(elemental_error);
2188 
2189  // Reset the minimum number of elements that need to be unrefined
2190  // to make it worthwhile
2191  max_keep_unrefined() = backup;
2192 
2193  // Now restore the error tolerances
2194  this->min_permitted_error() = current_min_error;
2195  this->max_permitted_error() = current_max_error;
2196 
2197  // Has the unrefinement actually changed anything?
2198  if (Nelement == this->nelement())
2199  {
2200  return 1;
2201  }
2202  else
2203  {
2204  return 0;
2205  }
2206  }
oomph::TreeBasedRefineableMeshBase::adapt
void adapt(const Vector< double > &elemental_error)
Definition: refineable_mesh.cc:307

References adapt(), e(), calibrate::error, oomph::RefineableMeshBase::max_keep_unrefined(), oomph::RefineableMeshBase::max_permitted_error(), oomph::RefineableMeshBase::min_permitted_error(), and oomph::Mesh::nelement().

Referenced by refine_base_mesh().

Member Data Documentation

◆ Forest_pt

TreeForest* oomph::TreeBasedRefineableMeshBase::Forest_pt
protected

Forest representation of the mesh.

Referenced by adapt_mesh(), oomph::ElasticRefineableQuarterPipeMesh< ELEMENT >::ElasticRefineableQuarterPipeMesh(), forest_pt(), oomph::TreeBasedRefineableMesh< ELEMENT >::p_refine_elements_if_required(), oomph::RefineableEighthSphereMesh< ELEMENT >::RefineableEighthSphereMesh(), oomph::RefineableFullCircleMesh< ELEMENT >::RefineableFullCircleMesh(), oomph::RefineableQuarterPipeMesh< ELEMENT >::RefineableQuarterPipeMesh(), oomph::RefineableQuarterTubeMesh< ELEMENT >::RefineableQuarterTubeMesh(), oomph::RefineableTubeMesh< ELEMENT >::RefineableTubeMesh(), oomph::RefineableLineMesh< ELEMENT >::setup_binary_tree_forest(), oomph::RefineableBrickMesh< ELEMENT >::setup_octree_forest(), oomph::RefineableQuadMesh< ELEMENT >::setup_quadtree_forest(), oomph::TreeBasedRefineableMesh< ELEMENT >::split_elements_if_required(), TreeBasedRefineableMeshBase(), oomph::RefineableQuarterPipeMesh< ELEMENT >::~RefineableQuarterPipeMesh(), and ~TreeBasedRefineableMeshBase().

◆ Max_p_refinement_level

unsigned oomph::TreeBasedRefineableMeshBase::Max_p_refinement_level
protected

Max. permissible p-refinement level (relative to base mesh)

Referenced by doc_adaptivity_targets(), max_p_refinement_level(), and TreeBasedRefineableMeshBase().

◆ Max_refinement_level

unsigned oomph::TreeBasedRefineableMeshBase::Max_refinement_level
protected

Max. permissible refinement level (relative to base mesh)

Referenced by doc_adaptivity_targets(), max_refinement_level(), and TreeBasedRefineableMeshBase().

◆ Min_p_refinement_level

unsigned oomph::TreeBasedRefineableMeshBase::Min_p_refinement_level
protected

Min. permissible p-refinement level (relative to base mesh)

Referenced by doc_adaptivity_targets(), min_p_refinement_level(), and TreeBasedRefineableMeshBase().

◆ Min_refinement_level

unsigned oomph::TreeBasedRefineableMeshBase::Min_refinement_level
protected

Min. permissible refinement level (relative to base mesh)

Referenced by doc_adaptivity_targets(), min_refinement_level(), and TreeBasedRefineableMeshBase().

◆ Uniform_refinement_level_when_pruned

unsigned oomph::TreeBasedRefineableMeshBase::Uniform_refinement_level_when_pruned
protected

Level to which the mesh was uniformly refined when it was pruned.

Referenced by TreeBasedRefineableMeshBase(), and uniform_refinement_level_when_pruned().

The documentation for this class was generated from the following files: