refineable_mesh.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #ifndef OOMPH_REFINEABLE_MESH_HEADER
27 #define OOMPH_REFINEABLE_MESH_HEADER
28 
29 #include "mesh.h"
30 #include "refineable_elements.h"
31 // Include the tree template to fill in the C++ split function
32 // Must be called after refineable_element.h
33 #include "tree.template.cc"
34 #include "error_estimator.h"
35 
36 namespace oomph
37 {
38  //=======================================================================
39  /// Base class for refineable meshes. Provides standardised interfaces
40  /// for the following standard mesh adaptation routines:
41  /// - Uniform mesh refinement.
42  /// - Adaptation, based on the elemental error estimates obtained
43  /// from an error estimator function which is accessed via a function
44  /// pointer.
45  ///
46  //=======================================================================
47  class RefineableMeshBase : public virtual Mesh
48  {
49  public:
50  bool adapt_flag()
51  {
52  return Adapt_flag;
53  }
54 
55  /// Constructor sets default values for refinement targets etc.
56  /// and initialises pointer to spatial error estimator to NULL.
58  {
59  // Initialise pointer to spatial error estimator
61 
62  // Targets
63  Min_permitted_error = 1.0e-5;
64  Max_permitted_error = 1.0e-3;
65 
66  // Actual errors
67  Min_error = 0.0;
68  Max_error = 0.0;
69 
70  // Do I adapt or not?
71  Adapt_flag = true;
72 
73  // Do I p-adapt or not?
74  P_adapt_flag = true;
75 
76  // Do I disable additional synchronisation of hanging nodes?
78 
79  // Where do I write the documatation of the refinement process?
80  Doc_info_pt = 0;
81 
82  // If only a few elements are scheduled for unrefinement, don't bother
83  // By default unrefine all
85 
86  // Number of elements where the refinement is over-ruled
88  };
89 
90 
91  /// Broken copy constructor
92  RefineableMeshBase(const RefineableMeshBase& dummy) = delete;
93 
94  /// Broken assignment operator
95  void operator=(const RefineableMeshBase&) = delete;
96 
97  /// Empty Destructor:
98  virtual ~RefineableMeshBase() {}
99 
100  /// Access fct for number of elements that were refined
101  unsigned nrefined()
102  {
103  return Nrefined;
104  }
105 
106  /// Access fct for number of elements that were unrefined
107  unsigned nunrefined()
108  {
109  return Nunrefined;
110  }
111 
112  /// Number of elements that would have liked to be refined further
113  /// but can't because they've reached the max. refinement level
115  {
116  return Nrefinement_overruled;
117  }
118 
119  /// Max. number of elements that we allow to remain unrefined
120  /// if no other mesh adaptation is required (to avoid
121  /// mesh-adaptations that would only unrefine a few elements
122  /// and then force a new solve -- this can't be worth our while!)
123  unsigned& max_keep_unrefined()
124  {
125  return Max_keep_unrefined;
126  }
127 
128  /// Doc the targets for mesh adaptation
129  virtual void doc_adaptivity_targets(std::ostream& outfile)
130  {
131  outfile << std::endl;
132  outfile << "Targets for mesh adaptation: " << std::endl;
133  outfile << "---------------------------- " << std::endl;
134  outfile << "Target for max. error: " << Max_permitted_error << std::endl;
135  outfile << "Target for min. error: " << Min_permitted_error << std::endl;
136  outfile << "Don't unrefine if less than " << Max_keep_unrefined
137  << " elements need unrefinement." << std::endl;
138  outfile << std::endl;
139  }
140 
141 
142  /// Access to spatial error estimator
144  {
146  }
147 
148  /// Access to spatial error estimator (const version
150  {
152  }
153 
154  /// Access fct for min. error (i.e. (try to) merge elements if
155  /// their error is smaller)
157  {
158  return Min_permitted_error;
159  }
160 
161  /// Access fct for max. error (i.e. split elements if their
162  /// error is larger)
164  {
165  return Max_permitted_error;
166  }
167 
168  /// Access fct for min. actual error in present solution (i.e. before
169  /// re-solve on adapted mesh)
170  double& min_error()
171  {
172  return Min_error;
173  }
174 
175  /// Access fct for max. actual error in present solution (i.e. before
176  /// re-solve on adapted mesh)
177  double& max_error()
178  {
179  return Max_error;
180  }
181 
182  /// Access fct for pointer to DocInfo
184  {
185  return Doc_info_pt;
186  }
187 
188  /// Enable adaptation
190  {
191  Adapt_flag = true;
192  }
193 
194  /// Disable adaptation
196  {
197  Adapt_flag = false;
198  }
199 
200  /// Enable adaptation
202  {
203  P_adapt_flag = true;
204  }
205 
206  /// Disable adaptation
208  {
209  P_adapt_flag = false;
210  }
211 
212  /// Enable additional synchronisation of hanging nodes
214  {
216  }
217 
218  /// Disable additional synchronisation of hanging nodes
220  {
222  }
223 
224  /// Return whether the mesh is to be adapted
226  {
227  return Adapt_flag;
228  }
229 
230  /// Return whether the mesh is to be adapted
232  {
233  return P_adapt_flag;
234  }
235 
236  /// Return whether additional synchronisation is enabled
238  {
240  }
241 
242  /// Access fct for DocInfo
244  {
245  return *Doc_info_pt;
246  }
247 
248  /// Adapt mesh: Refine elements whose error is lager than err_max
249  /// and (try to) unrefine those whose error is smaller than err_min
250  virtual void adapt(const Vector<double>& elemental_error) = 0;
251 
252  /// p-adapt mesh: Refine elements whose error is lager than err_max
253  /// and (try to) unrefine those whose error is smaller than err_min
254  virtual void p_adapt(const Vector<double>& elemental_error)
255  {
256  // Derived classes must implement this as required. Default throws an
257  // error.
258  std::ostringstream err_stream;
259  err_stream << "p_adapt() called in base class RefineableMeshBase."
260  << std::endl
261  << "This needs to be implemented in the derived class."
262  << std::endl;
263  throw OomphLibError(
264  err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
265  }
266 
267  /// Refine mesh uniformly and doc process
268  virtual void refine_uniformly(DocInfo& doc_info) = 0;
269 
270  /// Refine mesh uniformly
271  virtual void refine_uniformly()
272  {
274  doc_info.directory() = "";
277  }
278 
279  /// p-refine mesh uniformly and doc process
281  {
282  // Derived classes must implement this as required. Default throws an
283  // error.
284  std::ostringstream err_stream;
285  err_stream
286  << "p_refine_uniformly() called in base class RefineableMeshBase."
287  << std::endl
288  << "This needs to be implemented in the derived class." << std::endl;
289  throw OomphLibError(
290  err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
291  }
292 
293  /// p-refine mesh uniformly
294  virtual void p_refine_uniformly()
295  {
297  doc_info.directory() = "";
300  }
301 
302  /// Unrefine mesh uniformly: Return 0 for success,
303  /// 1 for failure (if unrefinement has reached the coarsest permitted
304  /// level)
305  virtual unsigned unrefine_uniformly() = 0;
306 
307  /// p-unrefine mesh uniformly
309  {
310  // Derived classes must implement this as required. Default throws an
311  // error.
312  std::ostringstream err_stream;
313  err_stream
314  << "p_unrefine_uniformly() called in base class RefineableMeshBase."
315  << std::endl
316  << "This needs to be implemented in the derived class." << std::endl;
317  throw OomphLibError(
318  err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
319  }
320 
321  protected:
322  /// Pointer to spatial error estimator
324 
325  /// Max. error (i.e. split elements if their error is larger)
327 
328  /// Min. error (i.e. (try to) merge elements if their error is smaller)
330 
331  /// Min.actual error
332  double Min_error;
333 
334  /// Max. actual error
335  double Max_error;
336 
337  /// Stats: Number of elements that were refined
338  unsigned Nrefined;
339 
340  /// Stats: Number of elements that were unrefined
341  unsigned Nunrefined;
342 
343  /// Flag that requests adaptation
345 
346  /// Flag that requests p-adaptation
348 
349  /// Flag that disables additional synchronisation of hanging nodes
351 
352  /// Pointer to DocInfo
354 
355  /// Max. number of elements that can remain unrefined
356  /// if no other mesh adaptation is required (to avoid
357  /// mesh-adaptations that would only unrefine a few elements
358  /// and then force a new solve -- this can't be worth our while!)
360 
361  /// Number of elements that would like to be refined further but
362  /// can't because they've reached the max. refinement level
364  };
365 
366 
367  /// /////////////////////////////////////////////////////////////////
368  /// /////////////////////////////////////////////////////////////////
369  /// /////////////////////////////////////////////////////////////////
370 
371 
372  //=======================================================================
373  /// Base class for tree-based refineable meshes.
374  //=======================================================================
376  {
377  public:
378  /// Constructor
380  {
381  // Max/min refinement levels
384 
385  // Max/min p-refinement levels
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
401  }
402 
403 
404  /// Broken copy constructor
406  delete;
407 
408  /// Broken assignment operator
409  void operator=(const TreeBasedRefineableMeshBase&) = delete;
410 
411  /// 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  }
421 
422  /// Adapt mesh: Refine elements whose error is lager than err_max
423  /// and (try to) unrefine those whose error is smaller than err_min
424  void adapt(const Vector<double>& elemental_error);
425 
426  /// p-adapt mesh: Refine elements whose error is lager than err_max
427  /// and (try to) unrefine those whose error is smaller than err_min
428  void p_adapt(const Vector<double>& elemental_error);
429 
430  /// Refine mesh uniformly and doc process
432 
433  /// Refine mesh uniformly
435  {
437  }
438 
439  /// p-refine mesh uniformly and doc process
441 
442  /// p-refine mesh uniformly
444  {
446  }
447 
448  /// Unrefine mesh uniformly: Return 0 for success,
449  /// 1 for failure (if unrefinement has reached the coarsest permitted
450  /// level)
451  unsigned unrefine_uniformly();
452 
453  /// p-unrefine mesh uniformly
455 
456  /// Set up the tree forest associated with the Mesh (if any)
457  virtual void setup_tree_forest() = 0;
458 
459  /// Return pointer to the Forest represenation of the mesh
461  {
462  return Forest_pt;
463  }
464 
465 
466  /// Doc the targets for mesh adaptation
467  void doc_adaptivity_targets(std::ostream& outfile)
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  }
484 
485  /// Access fct for max. permissible refinement level (relative to base mesh)
487  {
488  return Max_refinement_level;
489  }
490 
491  /// Access fct for min. permissible refinement level (relative to base mesh)
493  {
494  return Min_refinement_level;
495  }
496 
497  /// Access fct for max. permissible p-refinement level (relative to base
498  /// mesh)
500  {
501  return Max_p_refinement_level;
502  }
503 
504  /// Access fct for min. permissible p-refinement level (relative to base
505  /// mesh)
507  {
508  return Min_p_refinement_level;
509  }
510 
511  /// Perform the actual tree-based mesh adaptation,
512  /// documenting the progress in the directory specified in DocInfo object.
513  virtual void adapt_mesh(DocInfo& doc_info);
514 
515  /// Perform the actual tree-based mesh adaptation. A simple wrapper
516  /// to call the function without documentation.
517  virtual void adapt_mesh()
518  {
519  // Create a dummy doc_info object
521  doc_info.directory() = "";
523  // Call the other adapt mesh
525  }
526 
527  /// Perform the actual tree-based mesh p-adaptation,
528  /// documenting the progress in the directory specified in DocInfo object.
530 
531  /// Perform the actual tree-based mesh p-adaptation. A simple wrapper
532  /// to call the function without documentation.
534  {
535  // Create a dummy doc_info object
537  doc_info.directory() = "";
539  // Call the other adapt mesh
541  }
542 
543  /// Refine mesh by splitting the elements identified
544  /// by their numbers.
545  virtual void refine_selected_elements(
546  const Vector<unsigned>& elements_to_be_refined);
547 
548  /// Refine mesh by splitting the elements identified
549  /// by their pointers.
550  virtual void refine_selected_elements(
551  const Vector<RefineableElement*>& elements_to_be_refined);
552 
553  /// p-refine mesh by refining the elements identified
554  /// by their numbers.
556  const Vector<unsigned>& elements_to_be_refined);
557 
558  /// p-refine mesh by refining the elements identified
559  /// by their pointers.
561  const Vector<PRefineableElement*>& elements_to_be_refined_pt);
562 
563 
564  /// Refine base mesh to same degree as reference mesh (relative
565  /// to original unrefined mesh).
567  TreeBasedRefineableMeshBase* const& ref_mesh_pt);
568 
569  /// Refine base mesh to same degree as reference mesh minus one
570  /// level of refinement (relative to original unrefined mesh). Useful
571  /// function for multigrid solvers; allows the easy copy of a mesh
572  /// to the level of refinement just below the current one
574  TreeBasedRefineableMeshBase* const& ref_mesh_pt);
575 
576  /// Refine mesh once so that its topology etc becomes that of the
577  /// (finer!) reference mesh -- if possible! Useful for meshes in multigrid
578  /// hierarchies. If the meshes are too different and the conversion
579  /// cannot be performed, the code dies (provided PARANOID is enabled).
580  virtual void refine_as_in_reference_mesh(
581  TreeBasedRefineableMeshBase* const& ref_mesh_pt);
582 
583  /// Get max/min refinement levels in mesh
584  virtual void get_refinement_levels(unsigned& min_refinement_level,
585  unsigned& max_refinement_level);
586 
587  /// Extract the elements at a particular refinement level in
588  /// the refinement pattern - used in Mesh::redistribute or whatever it's
589  /// going to be called (RefineableMeshBase::reduce_halo_layers or something)
591  unsigned& refinement_level, Vector<RefineableElement*>& level_elements);
592 
593  /// Extract refinement pattern: Consider the hypothetical mesh
594  /// obtained by truncating the refinement of the current mesh to a given
595  /// level (where \c level=0 is the un-refined base mesh). To advance
596  /// to the next refinement level, we need to refine (split) the
597  /// \c to_be_refined[level].size() elements identified by the
598  /// element numbers contained in \c vector to_be_refined[level][...]
599  virtual void get_refinement_pattern(
600  Vector<Vector<unsigned>>& to_be_refined);
601 
602  /// Refine base mesh according to specified refinement pattern
603  void refine_base_mesh(Vector<Vector<unsigned>>& to_be_refined);
604 
605  /// Refine mesh according to refinement pattern in restart file
606  virtual void refine(std::ifstream& restart_file);
607 
608  /// Dump refinement pattern to allow for rebuild
609  virtual void dump_refinement(std::ostream& outfile);
610 
611  /// Read refinement pattern to allow for rebuild
612  virtual void read_refinement(std::ifstream& restart_file,
613  Vector<Vector<unsigned>>& to_be_refined);
614 
615 
616  /// Level to which the mesh was uniformly refined when it was pruned
617  /// (const version)
619  {
621  }
622 
623 
624  /// Level to which the mesh was uniformly refined when it was pruned
626  {
628  }
629 
630 #ifdef OOMPH_HAS_MPI
631 
632  /// Classify all halo and haloed information in the mesh (overloaded
633  /// version from Mesh base class. Calls that one first, then synchronises
634  /// hanging nodes)
636  const bool& report_stats)
637  {
638  // Call version in base class but don't bother to call
639  // resize_halo_nodes() -- we'll do it ourselves below
640  bool backup = Resize_halo_nodes_not_required;
644 
645  double t_start = 0.0;
647  {
648  t_start = TimingHelpers::timer();
649  }
650 
651  // Communicate positions of non-hanging nodes that depend on non-existent
652  // neighbours (e.g. h-refinement of neighbouring elements with different
653  // p-orders where the shared edge is on the outer edge of the halo layer)
655 
656  // Get number of continously interpolated variables
657  unsigned local_ncont_interpolated_values = 0;
658 
659  // Bypass if we have no elements and then get overall
660  // value by reduction
661  if (nelement() > 0)
662  {
663  local_ncont_interpolated_values =
664  dynamic_cast<RefineableElement*>(element_pt(0))
665  ->ncont_interpolated_values();
666  }
667  unsigned ncont_interpolated_values = 0;
668  MPI_Allreduce(&local_ncont_interpolated_values,
669  &ncont_interpolated_values,
670  1,
671  MPI_UNSIGNED,
672  MPI_MAX,
673  Comm_pt->mpi_comm());
674 
675  // Synchronise the hanging nodes
676  synchronise_hanging_nodes(ncont_interpolated_values);
677 
679  {
680  double t_end = TimingHelpers::timer();
681  oomph_info
682  << "Time for "
683  << "TreeBasedRefineableMeshBase::synchronise_hanging_nodes() "
684  << "incl. time for initial allreduce in "
685  << "TreeBasedRefineableMeshBase::classify_halo_and_haloed_nodes(): "
686  << t_end - t_start << std::endl;
687  }
688 
689  // Now do resize halo nodes after all (unless it's been
690  // declared to be unnecessary by somebody...)
692  {
694  }
695  }
696 
697 
698  /// Classify the halo and haloed nodes in the mesh.
699  void classify_halo_and_haloed_nodes(const bool& report_stats = false)
700  {
704  }
705 
706 #endif
707 
708  protected:
709 #ifdef OOMPH_HAS_MPI
710 
711  /// Synchronise the hanging nodes if the mesh is distributed
712  void synchronise_hanging_nodes(const unsigned& ncont_interpolated_values);
713 
714  /// Additional synchronisation of hanging nodes
715  /// Required for reconcilliation of hanging nodes on the outer edge of the
716  /// halo layer when using elements with nonuniformly spaced nodes
717  /// Must be implemented in templated derived class because ELEMENT template
718  /// parameter is required for the node-creation functions in the
719  /// Missing_masters_functions namespace
721  const unsigned& ncont_interpolated_values) = 0;
722 
723  /// Synchronise the positions of non-hanging nodes that depend on
724  /// non-existent neighbours (e.g. h-refinement of neighbouring elements
725  /// with different p-orders where the shared edge is on the outer edge of
726  /// the halo layer)
728 
729 #endif
730 
731  /// Split all the elements in the mesh if required. This template
732  /// free interface will be overloaded in RefineableMesh<ELEMENT> so that any
733  /// new elements that are created will be of the correct type.
734  virtual void split_elements_if_required() = 0;
735 
736  /// p-refine all the elements in the mesh if required. This template
737  /// free interface will be overloaded in RefineableMesh<ELEMENT> so that
738  /// any temporary copies of the element that are created will be of the
739  /// correct type.
740  virtual void p_refine_elements_if_required() = 0;
741 
742  /// Complete the hanging node scheme recursively
743  void complete_hanging_nodes(const int& ncont_interpolated_values);
744 
745 
746  /// Auxiliary routine for recursive hanging node completion
748  Vector<Node*>& master_nodes,
749  Vector<double>& hang_weights,
750  const int& ival);
751 
752  /// Level to which the mesh was uniformly refined when it was pruned
754 
755  /// Max. permissible refinement level (relative to base mesh)
757 
758  /// Min. permissible refinement level (relative to base mesh)
760 
761  /// Max. permissible p-refinement level (relative to base mesh)
763 
764  /// Min. permissible p-refinement level (relative to base mesh)
766 
767  /// Forest representation of the mesh
769 
770  private:
771 #ifdef OOMPH_HAS_MPI
772 
773  /// Helper struct to collate data required during
774  /// TreeBasedRefineableMeshBase::synchronise_hanging_nodes
776  {
780  double Weight;
783  // Store pointer to node and index of value in case translation attempt
784  // fails (see synchronise_hanging_nodes())
786  int icont;
787  };
788 
789 #endif
790  };
791 
792 
793  /// /////////////////////////////////////////////////////////////////
794  /// /////////////////////////////////////////////////////////////////
795  /// /////////////////////////////////////////////////////////////////
796 
797 
798  //=======================================================================
799  /// Templated base class for refineable meshes. The use of the template
800  /// parameter is required only for creating new elements during mesh
801  /// adaptation. This class overloaded the template-free inteface to
802  /// the function split_elements_if_required() to make use of the template
803  /// parameter.
804  /// All refineable meshes should inherit directly from
805  /// TreeBasedRefineableMesh<ELEMENT>
806  //=======================================================================
807  template<class ELEMENT>
809  {
810  private:
811  /// Split all the elements if required. Overload the template-free
812  /// interface so that any new elements that are created
813  /// will be of the correct type.
815  {
816  // Find the number of trees in the forest
817  unsigned n_tree = this->Forest_pt->ntree();
818  // Loop over all "active" elements in the forest and split them
819  // if required
820  for (unsigned long e = 0; e < n_tree; e++)
821  {
823  &Tree::split_if_required<ELEMENT>);
824  }
825  }
826 
827  /// p-refine all the elements if required. Overload the template-free
828  /// interface so that any temporary copies of the element that are created
829  /// will be of the correct type.
831  {
832  // BENFLAG: Make a non const pointer to the mesh so it can be passed
833  // (HACK)
834  Mesh* mesh_pt = this;
835  // Find the number of trees in the forest
836  unsigned n_tree = this->Forest_pt->ntree();
837  // Loop over all "active" elements in the forest and p-refine them
838  // if required
839  for (unsigned long e = 0; e < n_tree; e++)
840  {
842  &Tree::p_refine_if_required<ELEMENT>, mesh_pt);
843  }
844  }
845 
846  protected:
847 #ifdef OOMPH_HAS_MPI
848 
849  /// Additional setup of shared node scheme
850  /// This is Required for reconcilliation of hanging nodes acrross processor
851  /// boundaries when using elements with nonuniformly spaced nodes.
852  /// ELEMENT template parameter is required so that
853  /// MacroElementNodeUpdateNodes which are added as external halo master
854  /// nodes can be made fully functional
856  const unsigned& ncont_interpolated_values);
857 
858 #endif
859 
860  public:
861  /// Constructor, call the constructor of the base class
863 
864  /// Broken copy constructor
866 
867  /// Empty virtual destructor
869  };
870 
871 
872  /// /////////////////////////////////////////////////////////////////
873  /// /////////////////////////////////////////////////////////////////
874  /// /////////////////////////////////////////////////////////////////
875 
876 
877  //=======================================================================
878  /// Base class for refineable tet meshes
879  //=======================================================================
881  {
882  public:
883  /// Max element size allowed during adaptation
885  {
886  return Max_element_size;
887  }
888 
889  /// Min element size allowed during adaptation
891  {
892  return Min_element_size;
893  }
894 
895  /// Min edge ratio before remesh gets triggered
897  {
899  }
900 
901 
902  /// Doc the targets for mesh adaptation
903  void doc_adaptivity_targets(std::ostream& outfile)
904  {
905  outfile << std::endl;
906  outfile << "Targets for mesh adaptation: " << std::endl;
907  outfile << "---------------------------- " << std::endl;
908  outfile << "Target for max. error: " << Max_permitted_error << std::endl;
909  outfile << "Target for min. error: " << Min_permitted_error << std::endl;
910  outfile << "Target max edge ratio: " << Max_permitted_edge_ratio
911  << std::endl;
912  outfile << "Min. allowed element size: " << Min_element_size << std::endl;
913  outfile << "Max. allowed element size: " << Max_element_size << std::endl;
914  outfile << "Don't unrefine if less than " << Max_keep_unrefined
915  << " elements need unrefinement." << std::endl;
916  outfile << std::endl;
917  }
918 
919 
920  /// Compute target volume based on the elements' error and the
921  /// error target; return max edge ratio
922  double compute_volume_target(const Vector<double>& elem_error,
923  Vector<double>& target_volume)
924  {
925  double max_edge_ratio = 0.0;
926  unsigned count_unrefined = 0;
927  unsigned count_refined = 0;
928  this->Nrefinement_overruled = 0;
929 
930  unsigned nel = this->nelement();
931  for (unsigned e = 0; e < nel; e++)
932  {
933  // Get element
934  FiniteElement* el_pt = this->finite_element_pt(e);
935 
936  // Calculate the volume of the element
937  double volume = el_pt->size();
938 
939  // Find the vertex coordinates
940  // (vertices are enumerated first)
941  double vertex[4][3];
942  for (unsigned n = 0; n < 4; ++n)
943  {
944  for (unsigned i = 0; i < 3; ++i)
945  {
946  vertex[n][i] = el_pt->node_pt(n)->x(i);
947  }
948  }
949 
950  // Compute the radius of the circumsphere of the tetrahedron
951  // Algorithm stolen from tetgen for consistency
952  DenseDoubleMatrix A(3);
953  for (unsigned i = 0; i < 3; ++i)
954  {
955  A(0, i) = vertex[1][i] - vertex[0][i];
956  A(1, i) = vertex[2][i] - vertex[0][i];
957  A(2, i) = vertex[3][i] - vertex[0][i];
958  }
959 
960  Vector<double> rhs(3);
961  // Compute the right hand side vector b (3x1).
962  for (unsigned i = 0; i < 3; ++i)
963  {
964  rhs[i] = 0.0;
965  for (unsigned k = 0; k < 3; ++k)
966  {
967  rhs[i] += A(i, k) * A(i, k);
968  }
969  rhs[i] *= 0.5;
970  }
971 
972  // Solve the linear system, in which the rhs is over-written with
973  // the solution
974  A.solve(rhs);
975 
976  // Calculate the circum-radius
977  double circum_radius =
978  sqrt(rhs[0] * rhs[0] + rhs[1] * rhs[1] + rhs[2] * rhs[2]);
979 
980  // Now find the shortest edge length
981  Vector<double> edge(3);
982  double min_length = DBL_MAX;
983  for (unsigned start = 0; start < 4; ++start)
984  {
985  for (unsigned end = start + 1; end < 4; ++end)
986  {
987  for (unsigned i = 0; i < 3; ++i)
988  {
989  edge[i] = vertex[start][i] - vertex[end][i];
990  }
991  double length =
992  sqrt(edge[0] * edge[0] + edge[1] * edge[1] + edge[2] * edge[2]);
993  if (length < min_length)
994  {
995  min_length = length;
996  }
997  }
998  }
999 
1000  // Now calculate the minimum edge ratio for this element
1001  double local_max_edge_ratio = circum_radius / min_length;
1002  if (local_max_edge_ratio > max_edge_ratio)
1003  {
1004  max_edge_ratio = local_max_edge_ratio;
1005  }
1006 
1007  // Mimick refinement in tree-based procedure: Target volumes
1008  // for elements that exceed permitted error is 1/4 of their
1009  // current volume, corresponding to a uniform sub-division.
1010  if (elem_error[e] > this->max_permitted_error())
1011  {
1012  target_volume[e] = std::max(volume / 4.0, Min_element_size);
1013  if (target_volume[e] != Min_element_size)
1014  {
1015  count_refined++;
1016  }
1017  else
1018  {
1019  this->Nrefinement_overruled++;
1020  }
1021  }
1022  else if (elem_error[e] < this->min_permitted_error())
1023  {
1024  target_volume[e] = std::min(4.0 * volume, Max_element_size);
1025  if (target_volume[e] != Max_element_size)
1026  {
1027  count_unrefined++;
1028  }
1029  }
1030  else
1031  {
1032  // Leave it alone
1033  target_volume[e] = std::max(volume, Min_element_size);
1034  }
1035 
1036  } // End of loop over elements
1037 
1038  // Tell everybody
1039  this->Nrefined = count_refined;
1040  this->Nunrefined = count_unrefined;
1041 
1042  return max_edge_ratio;
1043  }
1044 
1045  /// Max permitted element size
1047 
1048  /// Min permitted element size
1050 
1051  /// Max edge ratio before remesh gets triggered
1053  };
1054 
1055 
1056 } // namespace oomph
1057 
1058 #endif
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1271
Information for documentation of results: Directory and file number to enable output in the form RESL...
void disable_doc()
Disable documentation.
std::string directory() const
Output directory.
void solve(DoubleVector &rhs)
Complete LU solve (replaces matrix by its LU decomposition and overwrites RHS with solution)....
Definition: matrices.cc:50
Base class for spatial error estimators.
A general Finite Element class.
Definition: elements.h:1313
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition: elements.cc:4290
Class that contains data for hanging nodes.
Definition: nodes.h:742
A general mesh class.
Definition: mesh.h:67
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
virtual void classify_halo_and_haloed_nodes(DocInfo &doc_info, const bool &report_stats)
Classify the halo and haloed nodes in the mesh. Virtual so it can be overloaded to perform additional...
Definition: mesh.cc:3262
OomphCommunicator * Comm_pt
Pointer to communicator – set to NULL if mesh is not distributed.
Definition: mesh.h:119
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
Definition: mesh.h:460
bool Resize_halo_nodes_not_required
Set this to true to suppress resizing of halo nodes (at your own risk!)
Definition: mesh.h:144
void resize_halo_nodes()
Helper function that resizes halo nodes to the same size as their non-halo counterparts if required....
Definition: mesh.cc:4426
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
Base class for refineable meshes. Provides standardised interfaces for the following standard mesh ad...
virtual void p_adapt(const Vector< double > &elemental_error)
p-adapt mesh: Refine elements whose error is lager than err_max and (try to) unrefine those whose err...
bool P_adapt_flag
Flag that requests p-adaptation.
unsigned Nrefined
Stats: Number of elements that were refined.
void enable_p_adaptation()
Enable adaptation.
unsigned nunrefined()
Access fct for number of elements that were unrefined.
unsigned & max_keep_unrefined()
Max. number of elements that we allow to remain unrefined if no other mesh adaptation is required (to...
DocInfo doc_info()
Access fct for DocInfo.
RefineableMeshBase()
Constructor sets default values for refinement targets etc. and initialises pointer to spatial error ...
bool Adapt_flag
Flag that requests adaptation.
void operator=(const RefineableMeshBase &)=delete
Broken assignment operator.
void enable_adaptation()
Enable adaptation.
DocInfo * Doc_info_pt
Pointer to DocInfo.
ErrorEstimator * spatial_error_estimator_pt() const
Access to spatial error estimator (const version.
virtual void p_refine_uniformly(DocInfo &doc_info)
p-refine mesh uniformly and doc process
unsigned Max_keep_unrefined
Max. number of elements that can remain unrefined if no other mesh adaptation is required (to avoid m...
virtual unsigned unrefine_uniformly()=0
Unrefine mesh uniformly: Return 0 for success, 1 for failure (if unrefinement has reached the coarses...
double & min_error()
Access fct for min. actual error in present solution (i.e. before re-solve on adapted mesh)
unsigned Nunrefined
Stats: Number of elements that were unrefined.
virtual void p_refine_uniformly()
p-refine mesh uniformly
void p_unrefine_uniformly(DocInfo &doc_info)
p-unrefine mesh uniformly
bool Additional_synchronisation_of_hanging_nodes_not_required
Flag that disables additional synchronisation of hanging nodes.
DocInfo *& doc_info_pt()
Access fct for pointer to DocInfo.
void disable_p_adaptation()
Disable adaptation.
void disable_adaptation()
Disable adaptation.
void enable_additional_synchronisation_of_hanging_nodes()
Enable additional synchronisation of hanging nodes.
bool is_additional_synchronisation_of_hanging_nodes_disabled() const
Return whether additional synchronisation is enabled.
double & max_error()
Access fct for max. actual error in present solution (i.e. before re-solve on adapted mesh)
bool is_adaptation_enabled() const
Return whether the mesh is to be adapted.
ErrorEstimator * Spatial_error_estimator_pt
Pointer to spatial error estimator.
bool is_p_adaptation_enabled() const
Return whether the mesh is to be adapted.
virtual void doc_adaptivity_targets(std::ostream &outfile)
Doc the targets for mesh adaptation.
double Max_error
Max. actual error.
virtual void refine_uniformly()
Refine mesh uniformly.
double Min_error
Min.actual error.
virtual void adapt(const Vector< double > &elemental_error)=0
Adapt mesh: Refine elements whose error is lager than err_max and (try to) unrefine those whose error...
unsigned nrefined()
Access fct for number of elements that were refined.
RefineableMeshBase(const RefineableMeshBase &dummy)=delete
Broken copy constructor.
double Min_permitted_error
Min. error (i.e. (try to) merge elements if their error is smaller)
virtual ~RefineableMeshBase()
Empty Destructor:
double & min_permitted_error()
Access fct for min. error (i.e. (try to) merge elements if their error is smaller)
unsigned Nrefinement_overruled
Number of elements that would like to be refined further but can't because they've reached the max....
void disable_additional_synchronisation_of_hanging_nodes()
Disable additional synchronisation of hanging nodes.
double & max_permitted_error()
Access fct for max. error (i.e. split elements if their error is larger)
unsigned & nrefinement_overruled()
Number of elements that would have liked to be refined further but can't because they've reached the ...
double Max_permitted_error
Max. error (i.e. split elements if their error is larger)
ErrorEstimator *& spatial_error_estimator_pt()
Access to spatial error estimator.
virtual void refine_uniformly(DocInfo &doc_info)=0
Refine mesh uniformly and doc process.
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
void doc_adaptivity_targets(std::ostream &outfile)
Doc the targets for mesh adaptation.
double compute_volume_target(const Vector< double > &elem_error, Vector< double > &target_volume)
Compute target volume based on the elements' error and the error target; return max edge ratio.
double & max_element_size()
Max element size allowed during adaptation.
double & min_element_size()
Min element size allowed during adaptation.
double Max_element_size
Max permitted element size.
double Min_element_size
Min permitted element size.
double Max_permitted_edge_ratio
Max edge ratio before remesh gets triggered.
double & max_permitted_edge_ratio()
Min edge ratio before remesh gets triggered.
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
void synchronise_hanging_nodes(const unsigned &ncont_interpolated_values)
Synchronise the hanging nodes if the mesh is distributed.
virtual void get_elements_at_refinement_level(unsigned &refinement_level, Vector< RefineableElement * > &level_elements)
Extract the elements at a particular refinement level in the refinement pattern - used in Mesh::redis...
void classify_halo_and_haloed_nodes(const bool &report_stats=false)
Classify the halo and haloed nodes in the mesh.
void refine_base_mesh(Vector< Vector< unsigned >> &to_be_refined)
Refine base mesh according to specified refinement pattern.
virtual void refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
Refine mesh by splitting the elements identified by their numbers.
unsigned & max_refinement_level()
Access fct for max. permissible refinement level (relative to base mesh)
unsigned Uniform_refinement_level_when_pruned
Level to which the mesh was uniformly refined when it was pruned.
unsigned uniform_refinement_level_when_pruned() const
Level to which the mesh was uniformly refined when it was pruned (const version)
virtual void dump_refinement(std::ostream &outfile)
Dump refinement pattern to allow for rebuild.
unsigned Max_p_refinement_level
Max. permissible p-refinement level (relative to base mesh)
TreeForest * forest_pt()
Return pointer to the Forest represenation of the mesh.
unsigned & min_p_refinement_level()
Access fct for min. permissible p-refinement level (relative to base mesh)
virtual void read_refinement(std::ifstream &restart_file, Vector< Vector< unsigned >> &to_be_refined)
Read refinement pattern to allow for rebuild.
unsigned unrefine_uniformly()
Unrefine mesh uniformly: Return 0 for success, 1 for failure (if unrefinement has reached the coarses...
virtual void refine(std::ifstream &restart_file)
Refine mesh according to refinement pattern in restart file.
virtual void additional_synchronise_hanging_nodes(const unsigned &ncont_interpolated_values)=0
Additional synchronisation of hanging nodes Required for reconcilliation of hanging nodes on the oute...
virtual void adapt_mesh()
Perform the actual tree-based mesh adaptation. A simple wrapper to call the function without document...
virtual void refine_base_mesh_as_in_reference_mesh(TreeBasedRefineableMeshBase *const &ref_mesh_pt)
Refine base mesh to same degree as reference mesh (relative to original unrefined mesh).
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.
virtual void refine_as_in_reference_mesh(TreeBasedRefineableMeshBase *const &ref_mesh_pt)
Refine mesh once so that its topology etc becomes that of the (finer!) reference mesh – if possible!...
unsigned & uniform_refinement_level_when_pruned()
Level to which the mesh was uniformly refined when it was pruned.
void operator=(const TreeBasedRefineableMeshBase &)=delete
Broken assignment operator.
void p_refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
p-refine mesh by refining the elements identified by their numbers.
virtual void split_elements_if_required()=0
Split all the elements in the mesh if required. This template free interface will be overloaded in Re...
void doc_adaptivity_targets(std::ostream &outfile)
Doc the targets for mesh adaptation.
TreeBasedRefineableMeshBase(const TreeBasedRefineableMeshBase &dummy)=delete
Broken copy constructor.
void refine_uniformly()
Refine mesh uniformly.
void p_unrefine_uniformly(DocInfo &doc_info)
p-unrefine mesh uniformly
virtual void get_refinement_levels(unsigned &min_refinement_level, unsigned &max_refinement_level)
Get max/min refinement levels in mesh.
unsigned & max_p_refinement_level()
Access fct for max. permissible p-refinement level (relative to base mesh)
virtual void p_refine_elements_if_required()=0
p-refine all the elements in the mesh if required. This template free interface will be overloaded in...
void complete_hanging_nodes(const int &ncont_interpolated_values)
Complete the hanging node scheme recursively.
virtual void setup_tree_forest()=0
Set up the tree forest associated with the Mesh (if any)
TreeForest * Forest_pt
Forest representation of the mesh.
virtual ~TreeBasedRefineableMeshBase()
Empty Destructor:
unsigned Min_p_refinement_level
Min. permissible p-refinement level (relative to base mesh)
void synchronise_nonhanging_nodes()
Synchronise the positions of non-hanging nodes that depend on non-existent neighbours (e....
void adapt(const Vector< double > &elemental_error)
Adapt mesh: Refine elements whose error is lager than err_max and (try to) unrefine those whose error...
unsigned Max_refinement_level
Max. permissible refinement level (relative to base mesh)
void p_adapt(const Vector< double > &elemental_error)
p-adapt mesh: Refine elements whose error is lager than err_max and (try to) unrefine those whose err...
void p_adapt_mesh()
Perform the actual tree-based mesh p-adaptation. A simple wrapper to call the function without docume...
virtual void get_refinement_pattern(Vector< Vector< unsigned >> &to_be_refined)
Extract refinement pattern: Consider the hypothetical mesh obtained by truncating the refinement of t...
virtual bool refine_base_mesh_as_in_reference_mesh_minus_one(TreeBasedRefineableMeshBase *const &ref_mesh_pt)
Refine base mesh to same degree as reference mesh minus one level of refinement (relative to original...
void p_refine_uniformly()
p-refine mesh uniformly
void classify_halo_and_haloed_nodes(DocInfo &doc_info, const bool &report_stats)
Classify all halo and haloed information in the mesh (overloaded version from Mesh base class....
unsigned & min_refinement_level()
Access fct for min. permissible refinement level (relative to base mesh)
unsigned Min_refinement_level
Min. permissible refinement level (relative to base mesh)
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
virtual ~TreeBasedRefineableMesh()
Empty virtual destructor.
void additional_synchronise_hanging_nodes(const unsigned &ncont_interpolated_values)
Additional setup of shared node scheme This is Required for reconcilliation of hanging nodes acrross ...
TreeBasedRefineableMesh(const TreeBasedRefineableMesh &dummy)=delete
Broken copy constructor.
void split_elements_if_required()
Split all the elements if required. Overload the template-free interface so that any new elements tha...
TreeBasedRefineableMesh()
Constructor, call the constructor of the base class.
void p_refine_elements_if_required()
p-refine all the elements if required. Overload the template-free interface so that any temporary cop...
A TreeForest consists of a collection of TreeRoots. Each member tree can have neighbours in various e...
Definition: tree.h:409
unsigned ntree()
Number of trees in forest.
Definition: tree.h:460
TreeRoot * tree_pt(const unsigned &i) const
Return pointer to i-th tree in forest.
Definition: tree.h:466
void traverse_leaves(Tree::VoidMemberFctPt member_function)
Traverse the tree and execute void Tree member function member_function() only at its leaves.
Definition: tree.cc:207
void start(const unsigned &i)
(Re-)start i-th timer
bool Doc_comprehensive_timings
Global boolean to switch on comprehensive timing – can probably be declared const false when developm...
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...
Helper struct to collate data required during TreeBasedRefineableMeshBase::synchronise_hanging_nodes.