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-2022 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
36namespace 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:
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:
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 {
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!)
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
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
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 {
489 }
490
491 /// Access fct for min. permissible refinement level (relative to base mesh)
493 {
495 }
496
497 /// Access fct for max. permissible p-refinement level (relative to base
498 /// mesh)
500 {
502 }
503
504 /// Access fct for min. permissible p-refinement level (relative to base
505 /// mesh)
507 {
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();
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.
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
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
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition: elements.cc:4290
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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...
unsigned & max_keep_unrefined()
Max. number of elements that we allow to remain unrefined if no other mesh adaptation is required (to...
bool P_adapt_flag
Flag that requests p-adaptation.
unsigned Nrefined
Stats: Number of elements that were refined.
double & min_permitted_error()
Access fct for min. error (i.e. (try to) merge elements if their error is smaller)
void enable_p_adaptation()
Enable adaptation.
unsigned nunrefined()
Access fct for number of elements that were unrefined.
DocInfo *& doc_info_pt()
Access fct for pointer to DocInfo.
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 & nrefinement_overruled()
Number of elements that would have liked to be refined further but can't because they've reached the ...
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...
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.
void disable_p_adaptation()
Disable adaptation.
void disable_adaptation()
Disable adaptation.
void enable_additional_synchronisation_of_hanging_nodes()
Enable additional synchronisation of hanging nodes.
ErrorEstimator *& spatial_error_estimator_pt()
Access to spatial error estimator.
bool is_additional_synchronisation_of_hanging_nodes_disabled() const
Return whether additional synchronisation is enabled.
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_error()
Access fct for min. actual error in present solution (i.e. before re-solve on adapted mesh)
double Min_permitted_error
Min. error (i.e. (try to) merge elements if their error is smaller)
virtual ~RefineableMeshBase()
Empty Destructor:
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_error()
Access fct for max. actual error in present solution (i.e. before re-solve on adapted mesh)
double & max_permitted_error()
Access fct for max. error (i.e. split elements if their error is larger)
double Max_permitted_error
Max. error (i.e. split elements if their error is larger)
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 & min_element_size()
Min element size allowed during adaptation.
double & max_element_size()
Max element size allowed during adaptation.
double & max_permitted_edge_ratio()
Min edge ratio before remesh gets triggered.
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.
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
TreeForest * forest_pt()
Return pointer to the Forest represenation of the mesh.
void synchronise_hanging_nodes(const unsigned &ncont_interpolated_values)
Synchronise the hanging nodes if the mesh is distributed.
virtual void read_refinement(std::ifstream &restart_file, Vector< Vector< unsigned > > &to_be_refined)
Read refinement pattern to allow for rebuild.
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.
virtual void refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
Refine mesh by splitting the elements identified by their numbers.
unsigned & max_p_refinement_level()
Access fct for max. permissible p-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)
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...
unsigned & min_refinement_level()
Access fct for min. permissible refinement level (relative to base mesh)
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).
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...
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.
unsigned & max_refinement_level()
Access fct for max. permissible refinement level (relative to base mesh)
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!...
void refine_base_mesh(Vector< Vector< unsigned > > &to_be_refined)
Refine base mesh according to specified refinement pattern.
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.
unsigned & min_p_refinement_level()
Access fct for min. permissible p-refinement level (relative to base mesh)
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.
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 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
Min. 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.
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
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.