refineable_elements.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//====================================================================
// Header file for classes that define refineable element objects
 
// Include guard to prevent multiple inclusions of the header
#ifndef OOMPH_REFINEABLE_ELEMENTS_HEADER
#define OOMPH_REFINEABLE_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include "elements.h"
#include "tree.h"
 
namespace oomph
{
  class Mesh;
 
  //=======================================================================
  //======================================================================
  class RefineableElement : public virtual FiniteElement
  {
  protected:
    Tree* Tree_pt;
 
    unsigned Refine_level;
 
    bool To_be_refined;
 
    bool Refinement_is_enabled;
 
    bool Sons_to_be_unrefined;
 
    long Number;
 
    static double Max_integrity_tolerance;
 
    static void check_value_id(const int& n_continuously_interpolated_values,
                               const int& value_id);
 
 
    void assemble_local_to_eulerian_jacobian(
      const DShape& dpsids, DenseMatrix<double>& jacobian) const;
 
    void assemble_local_to_eulerian_jacobian2(
      const DShape& d2psids, DenseMatrix<double>& jacobian2) const;
 
    void assemble_eulerian_base_vectors(
      const DShape& dpsids, DenseMatrix<double>& interpolated_G) const;
 
    double local_to_eulerian_mapping_diagonal(
      const DShape& dpsids,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& inverse_jacobian) const;
 
  private:
    std::map<Node*, int>* Local_hang_eqn;
 
    std::map<Node*, unsigned> Shape_controlling_node_lookup;
 
  protected:
    void assign_hanging_local_eqn_numbers(const bool& store_local_dof_pt);
 
    virtual void fill_in_jacobian_from_nodal_by_fd(
      Vector<double>& residuals, DenseMatrix<double>& jacobian);
 
  public:
    RefineableElement()
      : FiniteElement(),
        Tree_pt(0),
        Refine_level(0),
        To_be_refined(false),
        Refinement_is_enabled(true),
        Sons_to_be_unrefined(false),
        Number(-1),
        Local_hang_eqn(0)
    {
    }
 
    // (The body is now in the cc file to keep the xlC compiler happy under AIX)
    virtual ~RefineableElement();
 
    RefineableElement(const RefineableElement&) = delete;
 
    void operator=(const RefineableElement&) = delete;
 
    Tree* tree_pt()
    {
      return Tree_pt;
    }
 
    void set_tree_pt(Tree* my_tree_pt)
    {
      Tree_pt = my_tree_pt;
    }
 
    virtual unsigned required_nsons() const
    {
      return 0;
    }
 
 
    bool refinement_is_enabled()
    {
      return Refinement_is_enabled;
    }
 
    void disable_refinement()
    {
      Refinement_is_enabled = false;
    }
 
    void enable_refinement()
    {
      Refinement_is_enabled = true;
    }
 
    template<class ELEMENT>
    void split(Vector<ELEMENT*>& son_pt) const
    {
      // Increase refinement level
      int son_refine_level = Refine_level + 1;
 
      // How many sons are to be constructed
      unsigned n_sons = required_nsons();
      // Resize the son pointer
      son_pt.resize(n_sons);
 
      // Loop over the sons and construct
      for (unsigned i = 0; i < n_sons; i++)
      {
        son_pt[i] = new ELEMENT;
        // Set the refinement level of the newly constructed son.
        son_pt[i]->set_refinement_level(son_refine_level);
      }
    }
 
 
    inline int local_hang_eqn(Node* const& node_pt, const unsigned& i)
    {
#ifdef RANGE_CHECKING
      if (i > ncont_interpolated_values())
      {
        std::ostringstream error_message;
        error_message << "Range Error: Value " << i
                      << " is not in the range (0,"
                      << ncont_interpolated_values() - 1 << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      return Local_hang_eqn[i][node_pt];
    }
 
    virtual void build(Mesh*& mesh_pt,
                       Vector<Node*>& new_node_pt,
                       bool& was_already_built,
                       std::ofstream& new_nodes_file) = 0;
 
    void set_refinement_level(const int& refine_level)
    {
      Refine_level = refine_level;
    }
 
    unsigned refinement_level() const
    {
      return Refine_level;
    }
 
    void select_for_refinement()
    {
      To_be_refined = true;
    }
 
    void deselect_for_refinement()
    {
      To_be_refined = false;
    }
 
    void select_sons_for_unrefinement()
    {
      Sons_to_be_unrefined = true;
    }
 
    void deselect_sons_for_unrefinement()
    {
      Sons_to_be_unrefined = false;
    }
 
    bool to_be_refined()
    {
      return To_be_refined;
    }
 
    bool sons_to_be_unrefined()
    {
      return Sons_to_be_unrefined;
    }
 
    virtual void rebuild_from_sons(Mesh*& mesh_pt) = 0;
 
    virtual void unbuild()
    {
      // Get pointer to father element
      RefineableElement* father_pt = father_element_pt();
      // If there is no father, nothing to do
      if (father_pt == 0)
      {
        return;
      }
 
      // Loop over all the nodes
      unsigned n_node = this->nnode();
      for (unsigned n = 0; n < n_node; n++)
      {
        // If any node in this element is in the father, it can't be deleted
        if (father_pt->get_node_number(this->node_pt(n)) >= 0)
        {
          node_pt(n)->set_non_obsolete();
        }
      }
    }
 
    virtual void deactivate_element();
 
    virtual bool nodes_built()
    {
      return node_pt(0) != 0;
    }
 
    long number() const
    {
      return Number;
    }
 
    void set_number(const long& mynumber)
    {
      Number = mynumber;
    }
 
    virtual unsigned ncont_interpolated_values() const = 0;
 
    virtual void get_interpolated_values(const Vector<double>& s,
                                         Vector<double>& values)
    {
      get_interpolated_values(0, s, values);
    }
 
    virtual void get_interpolated_values(const unsigned& t,
                                         const Vector<double>& s,
                                         Vector<double>& values) = 0;
 
    virtual Node* interpolating_node_pt(const unsigned& n, const int& value_id)
 
    {
      return node_pt(n);
    }
 
    virtual double local_one_d_fraction_of_interpolating_node(
      const unsigned& n1d, const unsigned& i, const int& value_id)
    {
      return local_one_d_fraction_of_node(n1d, i);
    }
 
    virtual Node* get_interpolating_node_at_local_coordinate(
      const Vector<double>& s, const int& value_id)
 
    {
      return get_node_at_local_coordinate(s);
    }
 
 
    virtual unsigned ninterpolating_node(const int& value_id)
    {
      return nnode();
    }
 
    virtual unsigned ninterpolating_node_1d(const int& value_id)
    {
      return nnode_1d();
    }
 
    virtual void interpolating_basis(const Vector<double>& s,
                                     Shape& psi,
                                     const int& value_id) const
    {
      shape(s, psi);
    }
 
    virtual void check_integrity(double& max_error) = 0;
 
    static double& max_integrity_tolerance()
    {
      return Max_integrity_tolerance;
    }
 
    void identify_field_data_for_interactions(
      std::set<std::pair<Data*, unsigned>>& paired_field_data);
 
 
    inline void assign_nodal_local_eqn_numbers(const bool& store_local_dof_pt)
    {
      FiniteElement::assign_nodal_local_eqn_numbers(store_local_dof_pt);
      assign_hanging_local_eqn_numbers(store_local_dof_pt);
    }
 
    virtual RefineableElement* root_element_pt()
    {
      // If there is no tree -- the element is its own root
      if (Tree_pt == 0)
      {
        return this;
      }
      // Otherwise it's the tree's root object
      else
      {
        return Tree_pt->root_pt()->object_pt();
      }
    }
 
    virtual RefineableElement* father_element_pt() const
    {
      // If we have no tree, we have no father
      if (Tree_pt == 0)
      {
        return 0;
      }
      else
      {
        // Otherwise get the father of the tree
        Tree* father_pt = Tree_pt->father_pt();
        // If the tree has no father then return null, no father
        if (father_pt == 0)
        {
          return 0;
        }
        else
        {
          return father_pt->object_pt();
        }
      }
    }
 
    void get_father_at_refinement_level(
      unsigned& refinement_level, RefineableElement*& father_at_reflevel_pt)
    {
      // Get the father in the tree (it shouldn't try to get a null Tree...)
      Tree* father_pt = Tree_pt->father_pt();
      // Get the refineable element associated with this father
      RefineableElement* father_el_pt =
        dynamic_cast<RefineableElement*>(father_pt->object_pt());
      // Get the refinement level
      unsigned level = father_el_pt->refinement_level();
      // If the level matches the required one then return, if not call again
      if (level == refinement_level)
      {
        father_at_reflevel_pt = father_el_pt;
      }
      else
      {
        // Recursive call
        father_el_pt->get_father_at_refinement_level(refinement_level,
                                                     father_at_reflevel_pt);
      }
    }
 
    virtual void initial_setup(Tree* const& adopted_father_pt = 0,
                               const unsigned& initial_p_order = 0)
    {
    }
 
    virtual void pre_build(Mesh*& mesh_pt, Vector<Node*>& new_node_pt) {}
 
    virtual void further_build() {}
 
    virtual void setup_hanging_nodes(Vector<std::ofstream*>& output_stream) {}
 
    virtual void further_setup_hanging_nodes() {}
 
    void get_dresidual_dnodal_coordinates(
      RankThreeTensor<double>& dresidual_dnodal_coordinates);
 
 
    unsigned nshape_controlling_nodes()
    {
      return Shape_controlling_node_lookup.size();
    }
 
    std::map<Node*, unsigned> shape_controlling_node_lookup()
    {
      return Shape_controlling_node_lookup;
    }
  };
 
 
 
 
  //======================================================================
  //======================================================================
  class PRefineableElement : public virtual RefineableElement
  {
  protected:
    unsigned P_order;
 
    bool To_be_p_refined;
 
    bool P_refinement_is_enabled;
 
    bool To_be_p_unrefined;
 
  public:
    PRefineableElement()
      : RefineableElement(),
        P_order(2),
        To_be_p_refined(false),
        P_refinement_is_enabled(true),
        To_be_p_unrefined(false)
    {
    }
 
    virtual ~PRefineableElement() {}
 
    PRefineableElement(const PRefineableElement&) = delete;
 
    void operator=(const PRefineableElement&) = delete;
 
    bool p_refinement_is_enabled()
    {
      return P_refinement_is_enabled;
    }
 
    void disable_p_refinement()
    {
      P_refinement_is_enabled = false;
    }
 
    void enable_p_refinement()
    {
      P_refinement_is_enabled = true;
    }
 
    unsigned& p_order()
    {
      return P_order;
    }
 
    unsigned p_order() const
    {
      return P_order;
    }
 
    virtual unsigned initial_p_order() const = 0;
 
    void select_for_p_refinement()
    {
      To_be_p_refined = true;
    }
 
    void deselect_for_p_refinement()
    {
      To_be_p_refined = false;
    }
 
    void select_for_p_unrefinement()
    {
      To_be_p_unrefined = true;
    }
 
    void deselect_for_p_unrefinement()
    {
      To_be_p_unrefined = false;
    }
 
    bool to_be_p_refined()
    {
      return To_be_p_refined;
    }
 
    bool to_be_p_unrefined()
    {
      return To_be_p_unrefined;
    }
 
    virtual void p_refine(const int& inc,
                          Mesh* const& mesh_pt,
                          GeneralisedElement* const& clone_pt) = 0;
 
    // Overload the nodes_built function to check every node
    bool nodes_built()
    {
      // Must check that EVERY node exists
      unsigned n_node = this->nnode();
      for (unsigned n = 0; n < n_node; n++)
      {
        if (this->node_pt(n) == 0)
        {
          return false;
        }
      }
      // If we get here then all the nodes are built
      return true;
    }
  };
 
 
 
 
  //=======================================================================
  //======================================================================
  class NonRefineableElementWithHangingNodes : public virtual RefineableElement
  {
  public:
    void build(Mesh*& mesh_pt,
               Vector<Node*>& new_node_pt,
               bool& was_already_built,
               std::ofstream& new_nodes_file)
    {
      std::ostringstream error_message;
      error_message << "This function is broken as it's only needed/used \n"
                    << "during \"proper\" refinement\n";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
 
    void get_interpolated_values(const Vector<double>& s,
                                 Vector<double>& values)
    {
      std::ostringstream error_message;
      error_message << "This function is broken as it's only needed/used \n"
                    << "during \"proper\" refinement\n";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void get_interpolated_values(const unsigned& t,
                                         const Vector<double>& s,
                                         Vector<double>& values)
    {
      std::ostringstream error_message;
      error_message << "This function is broken as it's only needed/used \n"
                    << "during \"proper\" refinement\n";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
 
    void check_integrity(double& max_error)
    {
      std::ostringstream error_message;
      error_message << "This function is broken as it's only needed/used \n"
                    << "during \"proper\" refinement\n";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    void rebuild_from_sons(Mesh*& mesh_pt)
    {
      std::ostringstream error_message;
      error_message << "This function is broken as it's only needed/used \n"
                    << "during \"proper\" refinement\n";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  };
 
 
 
 
  //=======================================================================
  //======================================================================
  class RefineableSolidElement : public virtual RefineableElement,
                                 public virtual SolidFiniteElement
  {
  private:
    std::map<Node*, DenseMatrix<int>> Local_position_hang_eqn;
 
 
    void assign_solid_hanging_local_eqn_numbers(const bool& store_local_dof_pt);
 
 
  protected:
    bool Use_undeformed_macro_element_for_new_lagrangian_coords;
 
    void assemble_local_to_lagrangian_jacobian(
      const DShape& dpsids, DenseMatrix<double>& jacobian) const;
 
    void assemble_local_to_lagrangian_jacobian2(
      const DShape& d2psids, DenseMatrix<double>& jacobian2) const;
 
    double local_to_lagrangian_mapping_diagonal(
      const DShape& dpsids,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& inverse_jacobian) const;
 
  public:
    RefineableSolidElement()
      : RefineableElement(),
        SolidFiniteElement(),
        Use_undeformed_macro_element_for_new_lagrangian_coords(false)
    {
    }
 
    virtual ~RefineableSolidElement() {}
 
    void assign_solid_local_eqn_numbers(const bool& store_local_dof_pt)
    {
      SolidFiniteElement::assign_solid_local_eqn_numbers(store_local_dof_pt);
      assign_solid_hanging_local_eqn_numbers(store_local_dof_pt);
    }
 
    unsigned ngeom_data() const;
 
    Data* geom_data_pt(const unsigned& j);
 
    void identify_geometric_data(std::set<Data*>& geometric_data_pt);
 
    void fill_in_jacobian_from_solid_position_by_fd(
      Vector<double>& residuals, DenseMatrix<double>& jacobian);
 
    bool is_undeformed_macro_element_used_for_new_lagrangian_coords() const
    {
      return Use_undeformed_macro_element_for_new_lagrangian_coords;
    }
 
    void enable_use_of_undeformed_macro_element_for_new_lagrangian_coords()
    {
      Use_undeformed_macro_element_for_new_lagrangian_coords = true;
    }
 
    void disable_use_of_undeformed_macro_element_for_new_lagrangian_coords()
    {
      Use_undeformed_macro_element_for_new_lagrangian_coords = false;
    }
 
    DenseMatrix<int>& local_position_hang_eqn(Node* const& node_pt)
    {
      return Local_position_hang_eqn[node_pt];
    }
 
    virtual void further_build()
    {
      Use_undeformed_macro_element_for_new_lagrangian_coords =
        dynamic_cast<RefineableSolidElement*>(father_element_pt())
          ->is_undeformed_macro_element_used_for_new_lagrangian_coords();
 
      RefineableElement::further_build();
    }
  };
 
 
 
 
  //=======================================================================
  //======================================================================
  class NonRefineableSolidElementWithHangingNodes
    : public virtual NonRefineableElementWithHangingNodes,
      public virtual RefineableSolidElement
  {
  };
 
 
 
 
} // namespace oomph
 
#endif