refineable_mesh.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
#ifndef OOMPH_REFINEABLE_MESH_HEADER
#define OOMPH_REFINEABLE_MESH_HEADER
 
#include "mesh.h"
#include "refineable_elements.h"
// Include the tree template to fill in the C++ split function
// Must be called after refineable_element.h
#include "tree.template.cc"
#include "error_estimator.h"
 
namespace oomph
{
  //=======================================================================
  //=======================================================================
  class RefineableMeshBase : public virtual Mesh
  {
  public:
    bool adapt_flag()
    {
      return Adapt_flag;
    }
 
    RefineableMeshBase()
    {
      // Initialise pointer to spatial error estimator
      Spatial_error_estimator_pt = 0;
 
      // Targets
      Min_permitted_error = 1.0e-5;
      Max_permitted_error = 1.0e-3;
 
      // Actual errors
      Min_error = 0.0;
      Max_error = 0.0;
 
      // Do I adapt or not?
      Adapt_flag = true;
 
      // Do I p-adapt or not?
      P_adapt_flag = true;
 
      // Do I disable additional synchronisation of hanging nodes?
      Additional_synchronisation_of_hanging_nodes_not_required = false;
 
      // Where do I write the documatation of the refinement process?
      Doc_info_pt = 0;
 
      // If only a few elements are scheduled for unrefinement, don't bother
      // By default unrefine all
      Max_keep_unrefined = 0;
 
      // Number of elements where the refinement is over-ruled
      Nrefinement_overruled = 0;
    };
 
 
    RefineableMeshBase(const RefineableMeshBase& dummy) = delete;
 
    void operator=(const RefineableMeshBase&) = delete;
 
    virtual ~RefineableMeshBase() {}
 
    unsigned nrefined()
    {
      return Nrefined;
    }
 
    unsigned nunrefined()
    {
      return Nunrefined;
    }
 
    unsigned& nrefinement_overruled()
    {
      return Nrefinement_overruled;
    }
 
    unsigned& max_keep_unrefined()
    {
      return Max_keep_unrefined;
    }
 
    virtual void doc_adaptivity_targets(std::ostream& outfile)
    {
      outfile << std::endl;
      outfile << "Targets for mesh adaptation: " << std::endl;
      outfile << "---------------------------- " << std::endl;
      outfile << "Target for max. error: " << Max_permitted_error << std::endl;
      outfile << "Target for min. error: " << Min_permitted_error << std::endl;
      outfile << "Don't unrefine if less than " << Max_keep_unrefined
              << " elements need unrefinement." << std::endl;
      outfile << std::endl;
    }
 
 
    ErrorEstimator*& spatial_error_estimator_pt()
    {
      return Spatial_error_estimator_pt;
    }
 
    ErrorEstimator* spatial_error_estimator_pt() const
    {
      return Spatial_error_estimator_pt;
    }
 
    double& min_permitted_error()
    {
      return Min_permitted_error;
    }
 
    double& max_permitted_error()
    {
      return Max_permitted_error;
    }
 
    double& min_error()
    {
      return Min_error;
    }
 
    double& max_error()
    {
      return Max_error;
    }
 
    DocInfo*& doc_info_pt()
    {
      return Doc_info_pt;
    }
 
    void enable_adaptation()
    {
      Adapt_flag = true;
    }
 
    void disable_adaptation()
    {
      Adapt_flag = false;
    }
 
    void enable_p_adaptation()
    {
      P_adapt_flag = true;
    }
 
    void disable_p_adaptation()
    {
      P_adapt_flag = false;
    }
 
    void enable_additional_synchronisation_of_hanging_nodes()
    {
      Additional_synchronisation_of_hanging_nodes_not_required = false;
    }
 
    void disable_additional_synchronisation_of_hanging_nodes()
    {
      Additional_synchronisation_of_hanging_nodes_not_required = true;
    }
 
    bool is_adaptation_enabled() const
    {
      return Adapt_flag;
    }
 
    bool is_p_adaptation_enabled() const
    {
      return P_adapt_flag;
    }
 
    bool is_additional_synchronisation_of_hanging_nodes_disabled() const
    {
      return Additional_synchronisation_of_hanging_nodes_not_required;
    }
 
    DocInfo doc_info()
    {
      return *Doc_info_pt;
    }
 
    virtual void adapt(const Vector<double>& elemental_error) = 0;
 
    virtual void p_adapt(const Vector<double>& elemental_error)
    {
      // Derived classes must implement this as required. Default throws an
      // error.
      std::ostringstream err_stream;
      err_stream << "p_adapt() called in base class RefineableMeshBase."
                 << std::endl
                 << "This needs to be implemented in the derived class."
                 << std::endl;
      throw OomphLibError(
        err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void refine_uniformly(DocInfo& doc_info) = 0;
 
    virtual void refine_uniformly()
    {
      DocInfo doc_info;
      doc_info.directory() = "";
      doc_info.disable_doc();
      refine_uniformly(doc_info);
    }
 
    virtual void p_refine_uniformly(DocInfo& doc_info)
    {
      // Derived classes must implement this as required. Default throws an
      // error.
      std::ostringstream err_stream;
      err_stream
        << "p_refine_uniformly() called in base class RefineableMeshBase."
        << std::endl
        << "This needs to be implemented in the derived class." << std::endl;
      throw OomphLibError(
        err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void p_refine_uniformly()
    {
      DocInfo doc_info;
      doc_info.directory() = "";
      doc_info.disable_doc();
      p_refine_uniformly(doc_info);
    }
 
    virtual unsigned unrefine_uniformly() = 0;
 
    void p_unrefine_uniformly(DocInfo& doc_info)
    {
      // Derived classes must implement this as required. Default throws an
      // error.
      std::ostringstream err_stream;
      err_stream
        << "p_unrefine_uniformly() called in base class RefineableMeshBase."
        << std::endl
        << "This needs to be implemented in the derived class." << std::endl;
      throw OomphLibError(
        err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
  protected:
    ErrorEstimator* Spatial_error_estimator_pt;
 
    double Max_permitted_error;
 
    double Min_permitted_error;
 
    double Min_error;
 
    double Max_error;
 
    unsigned Nrefined;
 
    unsigned Nunrefined;
 
    bool Adapt_flag;
 
    bool P_adapt_flag;
 
    bool Additional_synchronisation_of_hanging_nodes_not_required;
 
    DocInfo* Doc_info_pt;
 
    unsigned Max_keep_unrefined;
 
    unsigned Nrefinement_overruled;
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  class TreeBasedRefineableMeshBase : public virtual RefineableMeshBase
  {
  public:
    TreeBasedRefineableMeshBase()
    {
      // Max/min refinement levels
      Max_refinement_level = 5;
      Min_refinement_level = 0;
 
      // Max/min p-refinement levels
      Max_p_refinement_level = 7;
      Min_p_refinement_level = 2;
 
      // Stats
      Nrefined = 0;
      Nunrefined = 0;
 
      // Where do I write the documatation of the refinement process?
      Doc_info_pt = 0;
 
      // Initialise the forest pointer to NULL
      Forest_pt = 0;
 
      // Mesh hasn't been pruned yet
      Uniform_refinement_level_when_pruned = 0;
    }
 
 
    TreeBasedRefineableMeshBase(const TreeBasedRefineableMeshBase& dummy) =
      delete;
 
    void operator=(const TreeBasedRefineableMeshBase&) = delete;
 
    virtual ~TreeBasedRefineableMeshBase()
    {
      // Kill the forest if there is one
      if (Forest_pt != 0)
      {
        delete Forest_pt;
        Forest_pt = 0;
      }
    }
 
    void adapt(const Vector<double>& elemental_error);
 
    void p_adapt(const Vector<double>& elemental_error);
 
    void refine_uniformly(DocInfo& doc_info);
 
    void refine_uniformly()
    {
      RefineableMeshBase::refine_uniformly();
    }
 
    void p_refine_uniformly(DocInfo& doc_info);
 
    void p_refine_uniformly()
    {
      RefineableMeshBase::p_refine_uniformly();
    }
 
    unsigned unrefine_uniformly();
 
    void p_unrefine_uniformly(DocInfo& doc_info);
 
    virtual void setup_tree_forest() = 0;
 
    TreeForest* forest_pt()
    {
      return Forest_pt;
    }
 
 
    void doc_adaptivity_targets(std::ostream& outfile)
    {
      outfile << std::endl;
      outfile << "Targets for mesh adaptation: " << std::endl;
      outfile << "---------------------------- " << std::endl;
      outfile << "Target for max. error: " << Max_permitted_error << std::endl;
      outfile << "Target for min. error: " << Min_permitted_error << std::endl;
      outfile << "Min. refinement level: " << Min_refinement_level << std::endl;
      outfile << "Max. refinement level: " << Max_refinement_level << std::endl;
      outfile << "Min. p-refinement level: " << Min_p_refinement_level
              << std::endl;
      outfile << "Max. p-refinement level: " << Max_p_refinement_level
              << std::endl;
      outfile << "Don't unrefine if less than " << Max_keep_unrefined
              << " elements need unrefinement." << std::endl;
      outfile << std::endl;
    }
 
    unsigned& max_refinement_level()
    {
      return Max_refinement_level;
    }
 
    unsigned& min_refinement_level()
    {
      return Min_refinement_level;
    }
 
    unsigned& max_p_refinement_level()
    {
      return Max_p_refinement_level;
    }
 
    unsigned& min_p_refinement_level()
    {
      return Min_p_refinement_level;
    }
 
    virtual void adapt_mesh(DocInfo& doc_info);
 
    virtual void adapt_mesh()
    {
      // Create a dummy doc_info object
      DocInfo doc_info;
      doc_info.directory() = "";
      doc_info.disable_doc();
      // Call the other adapt mesh
      adapt_mesh(doc_info);
    }
 
    void p_adapt_mesh(DocInfo& doc_info);
 
    void p_adapt_mesh()
    {
      // Create a dummy doc_info object
      DocInfo doc_info;
      doc_info.directory() = "";
      doc_info.disable_doc();
      // Call the other adapt mesh
      p_adapt_mesh(doc_info);
    }
 
    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);
 
    void p_refine_selected_elements(
      const Vector<PRefineableElement*>& elements_to_be_refined_pt);
 
 
    virtual void refine_base_mesh_as_in_reference_mesh(
      TreeBasedRefineableMeshBase* const& ref_mesh_pt);
 
    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);
 
    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);
 
    virtual void refine(std::ifstream& restart_file);
 
    virtual void dump_refinement(std::ostream& outfile);
 
    virtual void read_refinement(std::ifstream& restart_file,
                                 Vector<Vector<unsigned>>& to_be_refined);
 
 
    unsigned uniform_refinement_level_when_pruned() const
    {
      return Uniform_refinement_level_when_pruned;
    }
 
 
    unsigned& uniform_refinement_level_when_pruned()
    {
      return Uniform_refinement_level_when_pruned;
    }
 
#ifdef OOMPH_HAS_MPI
 
    void classify_halo_and_haloed_nodes(DocInfo& doc_info,
                                        const bool& report_stats)
    {
      // Call version in base class but don't bother to call
      // resize_halo_nodes() -- we'll do it ourselves below
      bool backup = Resize_halo_nodes_not_required;
      Resize_halo_nodes_not_required = false;
      Mesh::classify_halo_and_haloed_nodes(doc_info, report_stats);
      Resize_halo_nodes_not_required = backup;
 
      double t_start = 0.0;
      if (Global_timings::Doc_comprehensive_timings)
      {
        t_start = TimingHelpers::timer();
      }
 
      // Communicate positions of non-hanging nodes that depend on non-existent
      // neighbours (e.g. h-refinement of neighbouring elements with different
      // p-orders where the shared edge is on the outer edge of the halo layer)
      synchronise_nonhanging_nodes();
 
      // Get number of continously interpolated variables
      unsigned local_ncont_interpolated_values = 0;
 
      // Bypass if we have no elements and then get overall
      // value by reduction
      if (nelement() > 0)
      {
        local_ncont_interpolated_values =
          dynamic_cast<RefineableElement*>(element_pt(0))
            ->ncont_interpolated_values();
      }
      unsigned ncont_interpolated_values = 0;
      MPI_Allreduce(&local_ncont_interpolated_values,
                    &ncont_interpolated_values,
                    1,
                    MPI_UNSIGNED,
                    MPI_MAX,
                    Comm_pt->mpi_comm());
 
      // Synchronise the hanging nodes
      synchronise_hanging_nodes(ncont_interpolated_values);
 
      if (Global_timings::Doc_comprehensive_timings)
      {
        double t_end = TimingHelpers::timer();
        oomph_info
          << "Time for "
          << "TreeBasedRefineableMeshBase::synchronise_hanging_nodes() "
          << "incl. time for initial allreduce in "
          << "TreeBasedRefineableMeshBase::classify_halo_and_haloed_nodes(): "
          << t_end - t_start << std::endl;
      }
 
      // Now do resize halo nodes after all (unless it's been
      // declared to be unnecessary by somebody...)
      if (!Resize_halo_nodes_not_required)
      {
        resize_halo_nodes();
      }
    }
 
 
    void classify_halo_and_haloed_nodes(const bool& report_stats = false)
    {
      DocInfo doc_info;
      doc_info.disable_doc();
      classify_halo_and_haloed_nodes(doc_info, report_stats);
    }
 
#endif
 
  protected:
#ifdef OOMPH_HAS_MPI
 
    void synchronise_hanging_nodes(const unsigned& ncont_interpolated_values);
 
    virtual void additional_synchronise_hanging_nodes(
      const unsigned& ncont_interpolated_values) = 0;
 
    void synchronise_nonhanging_nodes();
 
#endif
 
    virtual void split_elements_if_required() = 0;
 
    virtual void p_refine_elements_if_required() = 0;
 
    void complete_hanging_nodes(const int& ncont_interpolated_values);
 
 
    void complete_hanging_nodes_recursively(Node*& nod_pt,
                                            Vector<Node*>& master_nodes,
                                            Vector<double>& hang_weights,
                                            const int& ival);
 
    unsigned Uniform_refinement_level_when_pruned;
 
    unsigned Max_refinement_level;
 
    unsigned Min_refinement_level;
 
    unsigned Max_p_refinement_level;
 
    unsigned Min_p_refinement_level;
 
    TreeForest* Forest_pt;
 
  private:
#ifdef OOMPH_HAS_MPI
 
    struct HangHelperStruct
    {
      unsigned Sending_processor;
      unsigned Shared_node_id_on_sending_processor;
      unsigned Shared_node_proc;
      double Weight;
      HangInfo* Hang_pt;
      unsigned Master_node_index;
      // Store pointer to node and index of value in case translation attempt
      // fails (see synchronise_hanging_nodes())
      Node* Node_pt;
      int icont;
    };
 
#endif
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class TreeBasedRefineableMesh : public virtual TreeBasedRefineableMeshBase
  {
  private:
    void split_elements_if_required()
    {
      // Find the number of trees in the forest
      unsigned n_tree = this->Forest_pt->ntree();
      // Loop over all "active" elements in the forest and split them
      // if required
      for (unsigned long e = 0; e < n_tree; e++)
      {
        this->Forest_pt->tree_pt(e)->traverse_leaves(
          &Tree::split_if_required<ELEMENT>);
      }
    }
 
    void p_refine_elements_if_required()
    {
      // BENFLAG: Make a non const pointer to the mesh so it can be passed
      // (HACK)
      Mesh* mesh_pt = this;
      // Find the number of trees in the forest
      unsigned n_tree = this->Forest_pt->ntree();
      // Loop over all "active" elements in the forest and p-refine them
      // if required
      for (unsigned long e = 0; e < n_tree; e++)
      {
        this->Forest_pt->tree_pt(e)->traverse_leaves(
          &Tree::p_refine_if_required<ELEMENT>, mesh_pt);
      }
    }
 
  protected:
#ifdef OOMPH_HAS_MPI
 
    void additional_synchronise_hanging_nodes(
      const unsigned& ncont_interpolated_values);
 
#endif
 
  public:
    TreeBasedRefineableMesh() : TreeBasedRefineableMeshBase() {}
 
    TreeBasedRefineableMesh(const TreeBasedRefineableMesh& dummy) = delete;
 
    virtual ~TreeBasedRefineableMesh() {}
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  class RefineableTetMeshBase : public virtual RefineableMeshBase
  {
  public:
    double& max_element_size()
    {
      return Max_element_size;
    }
 
    double& min_element_size()
    {
      return Min_element_size;
    }
 
    double& max_permitted_edge_ratio()
    {
      return Max_permitted_edge_ratio;
    }
 
 
    void doc_adaptivity_targets(std::ostream& outfile)
    {
      outfile << std::endl;
      outfile << "Targets for mesh adaptation: " << std::endl;
      outfile << "---------------------------- " << std::endl;
      outfile << "Target for max. error: " << Max_permitted_error << std::endl;
      outfile << "Target for min. error: " << Min_permitted_error << std::endl;
      outfile << "Target max edge ratio: " << Max_permitted_edge_ratio
              << std::endl;
      outfile << "Min. allowed element size: " << Min_element_size << std::endl;
      outfile << "Max. allowed element size: " << Max_element_size << std::endl;
      outfile << "Don't unrefine if less than " << Max_keep_unrefined
              << " elements need unrefinement." << std::endl;
      outfile << std::endl;
    }
 
 
    double compute_volume_target(const Vector<double>& elem_error,
                                 Vector<double>& target_volume)
    {
      double max_edge_ratio = 0.0;
      unsigned count_unrefined = 0;
      unsigned count_refined = 0;
      this->Nrefinement_overruled = 0;
 
      unsigned nel = this->nelement();
      for (unsigned e = 0; e < nel; e++)
      {
        // Get element
        FiniteElement* el_pt = this->finite_element_pt(e);
 
        // Calculate the volume of the element
        double volume = el_pt->size();
 
        // Find the vertex coordinates
        // (vertices are enumerated first)
        double vertex[4][3];
        for (unsigned n = 0; n < 4; ++n)
        {
          for (unsigned i = 0; i < 3; ++i)
          {
            vertex[n][i] = el_pt->node_pt(n)->x(i);
          }
        }
 
        // Compute the radius of the circumsphere of the tetrahedron
        // Algorithm stolen from tetgen for consistency
        DenseDoubleMatrix A(3);
        for (unsigned i = 0; i < 3; ++i)
        {
          A(0, i) = vertex[1][i] - vertex[0][i];
          A(1, i) = vertex[2][i] - vertex[0][i];
          A(2, i) = vertex[3][i] - vertex[0][i];
        }
 
        Vector<double> rhs(3);
        // Compute the right hand side vector b (3x1).
        for (unsigned i = 0; i < 3; ++i)
        {
          rhs[i] = 0.0;
          for (unsigned k = 0; k < 3; ++k)
          {
            rhs[i] += A(i, k) * A(i, k);
          }
          rhs[i] *= 0.5;
        }
 
        // Solve the linear system, in which the rhs is over-written with
        // the solution
        A.solve(rhs);
 
        // Calculate the circum-radius
        double circum_radius =
          sqrt(rhs[0] * rhs[0] + rhs[1] * rhs[1] + rhs[2] * rhs[2]);
 
        // Now find the shortest edge length
        Vector<double> edge(3);
        double min_length = DBL_MAX;
        for (unsigned start = 0; start < 4; ++start)
        {
          for (unsigned end = start + 1; end < 4; ++end)
          {
            for (unsigned i = 0; i < 3; ++i)
            {
              edge[i] = vertex[start][i] - vertex[end][i];
            }
            double length =
              sqrt(edge[0] * edge[0] + edge[1] * edge[1] + edge[2] * edge[2]);
            if (length < min_length)
            {
              min_length = length;
            }
          }
        }
 
        // Now calculate the minimum edge ratio for this element
        double local_max_edge_ratio = circum_radius / min_length;
        if (local_max_edge_ratio > max_edge_ratio)
        {
          max_edge_ratio = local_max_edge_ratio;
        }
 
        // Mimick refinement in tree-based procedure: Target volumes
        // for elements that exceed permitted error is 1/4 of their
        // current volume, corresponding to a uniform sub-division.
        if (elem_error[e] > this->max_permitted_error())
        {
          target_volume[e] = std::max(volume / 4.0, Min_element_size);
          if (target_volume[e] != Min_element_size)
          {
            count_refined++;
          }
          else
          {
            this->Nrefinement_overruled++;
          }
        }
        else if (elem_error[e] < this->min_permitted_error())
        {
          target_volume[e] = std::min(4.0 * volume, Max_element_size);
          if (target_volume[e] != Max_element_size)
          {
            count_unrefined++;
          }
        }
        else
        {
          // Leave it alone
          target_volume[e] = std::max(volume, Min_element_size);
        }
 
      } // End of loop over elements
 
      // Tell everybody
      this->Nrefined = count_refined;
      this->Nunrefined = count_unrefined;
 
      return max_edge_ratio;
    }
 
    double Max_element_size;
 
    double Min_element_size;
 
    double Max_permitted_edge_ratio;
  };
 
 
} // namespace oomph
 
#endif