spines.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 spine nodes, elements and meshes
 
// Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_SPINES_HEADER
#define OOMPH_SPINES_HEADER
 
 
#include <string>
 
// oomph-lib headers
#include "nodes.h"
#include "elements.h"
#include "mesh.h"
#include "geom_objects.h"
#include "element_with_moving_nodes.h"
 
namespace oomph
{
  //=================================================================
  //=================================================================
  class Spine
  {
  public:
    Spine()
    {
      Geom_data_pt.resize(1);
      // Create Data for height. By default it's free
      Geom_data_pt[0] = new Data(1);
    }
 
    Spine(const double& height)
    {
      Geom_data_pt.resize(1);
      // Create Data for height. By default it's free
      Geom_data_pt[0] = new Data(1);
      // Set value
      Geom_data_pt[0]->set_value(0, height);
    }
 
    Spine(const double& height, const Vector<Data*>& geom_data_pt)
    {
      // Find the number of geometric data passed
      const unsigned n_geom_data = geom_data_pt.size();
      // Now allocate enough storage for the spine height and additional
      // geometric data
      Geom_data_pt.resize(n_geom_data + 1);
 
      // Create Data for height. By default it's free
      Geom_data_pt[0] = new Data(1);
      // Set value
      Geom_data_pt[0]->set_value(0, height);
      // Add the additional geometric data
      for (unsigned i = 0; i < n_geom_data; i++)
      {
        Geom_data_pt[i + 1] = geom_data_pt[i];
      }
    }
 
    Spine(const double& height,
          const Vector<Data*>& geom_data_pt,
          const Vector<GeomObject*>& geom_object_pt)
      : Geom_object_pt(geom_object_pt)
    {
      // Find the number of geometric data passed
      const unsigned n_geom_data = geom_data_pt.size();
      // Now allocate enough storage for the spine height and additional
      // geometric data
      Geom_data_pt.resize(n_geom_data + 1);
 
      // Create Data for height. By default it's free
      Geom_data_pt[0] = new Data(1);
      // Set value
      Geom_data_pt[0]->set_value(0, height);
      // Add the additional geometric data
      for (unsigned i = 0; i < n_geom_data; i++)
      {
        Geom_data_pt[i + 1] = geom_data_pt[i];
      }
    }
 
 
    ~Spine()
    {
      // Kill spine height
      delete Geom_data_pt[0];
    }
 
 
    double& height()
    {
      return *(Geom_data_pt[0]->value_pt(0));
    }
 
    Data*& spine_height_pt()
    {
      return Geom_data_pt[0];
    }
 
    Data* spine_height_pt() const
    {
      return Geom_data_pt[0];
    }
 
 
    unsigned ngeom_data()
    {
      return Geom_data_pt.size();
    }
 
    void set_geom_data_pt(const Vector<Data*>& geom_data_pt)
    {
      unsigned n_geom_data = geom_data_pt.size();
      Geom_data_pt.resize(n_geom_data + 1);
      for (unsigned i = 1; i < n_geom_data; i++)
      {
        Geom_data_pt[i + 1] = geom_data_pt[i];
      }
    }
 
    void add_geom_data_pt(Data* geom_data_pt)
    {
      Geom_data_pt.push_back(geom_data_pt);
    }
 
    Data*& geom_data_pt(const unsigned& i)
    {
      return Geom_data_pt[i];
    }
 
    Data* geom_data_pt(const unsigned& i) const
    {
      return Geom_data_pt[i];
    }
 
    Vector<Data*>& vector_geom_data_pt()
    {
      return Geom_data_pt;
    }
 
    unsigned ngeom_object()
    {
      return Geom_object_pt.size();
    }
 
    void set_geom_object_pt(const Vector<GeomObject*>& geom_object_pt)
    {
      unsigned n_geom_object = geom_object_pt.size();
      Geom_object_pt.resize(n_geom_object);
      for (unsigned i = 0; i < n_geom_object; i++)
      {
        Geom_object_pt[i] = geom_object_pt[i];
      }
    }
 
    void add_geom_object_pt(GeomObject* geom_object_pt)
    {
      Geom_object_pt.push_back(geom_object_pt);
    }
 
    GeomObject*& geom_object_pt(const unsigned& i)
    {
      return Geom_object_pt[i];
    }
 
    GeomObject* geom_object_pt(const unsigned& i) const
    {
      return Geom_object_pt[i];
    }
 
    Vector<GeomObject*>& vector_geom_object_pt()
    {
      return Geom_object_pt;
    }
 
    unsigned ngeom_parameter()
    {
      return Geom_parameter.size();
    }
 
    void set_geom_parameter(const Vector<double>& geom_parameter)
    {
      Geom_parameter = geom_parameter;
    }
 
    void add_geom_parameter(const double& geom_parameter)
    {
      Geom_parameter.push_back(geom_parameter);
    }
 
    double& geom_parameter(const unsigned& i)
    {
      return Geom_parameter[i];
    }
 
    const double& geom_parameter(const unsigned& i) const
    {
      return Geom_parameter[i];
    }
 
 
  private:
    // Data* Spine_height_pt;
 
    Vector<Data*> Geom_data_pt;
 
    Vector<GeomObject*> Geom_object_pt;
 
    Vector<double> Geom_parameter;
  };
 
 
  // Forward declaration
  class SpineMesh;
 
 
  //=====================================================================
  //=====================================================================
  class SpineNode : public Node
  {
  private:
    Spine* Spine_pt;
 
    double Fraction;
 
    SpineMesh* Spine_mesh_pt;
 
    unsigned Node_update_fct_id;
 
 
  public:
    SpineNode(const unsigned& n_dim,
              const unsigned& n_position_type,
              const unsigned& initial_nvalue)
      : Node(n_dim, n_position_type, initial_nvalue),
        Spine_pt(0),
        Fraction(0),
        Spine_mesh_pt(0),
        Node_update_fct_id(0)
    {
    }
 
    SpineNode(TimeStepper* const& time_stepper_pt,
              const unsigned& n_dim,
              const unsigned& n_position_type,
              const unsigned& initial_nvalue)
      : Node(time_stepper_pt, n_dim, n_position_type, initial_nvalue),
        Spine_pt(0),
        Fraction(0),
        Spine_mesh_pt(0),
        Node_update_fct_id(0)
    {
    }
 
    Spine*& spine_pt()
    {
      return Spine_pt;
    }
 
    double& fraction()
    {
      return Fraction;
    }
 
    unsigned& node_update_fct_id()
    {
      return Node_update_fct_id;
    }
 
    SpineMesh*& spine_mesh_pt()
    {
      return Spine_mesh_pt;
    }
 
    double& h()
    {
      return Spine_pt->height();
    }
 
    void node_update(const bool& update_all_time_levels_for_new_node = false);
 
    unsigned ngeom_data() const
    {
      if (Spine_pt)
      {
        return Spine_pt->ngeom_data();
      }
      else
      {
        return 0;
      }
    }
 
    unsigned ngeom_object() const
    {
      if (Spine_pt)
      {
        return Spine_pt->ngeom_object();
      }
      else
      {
        return 0;
      }
    }
 
    Data** all_geom_data_pt()
    {
      return &(Spine_pt->geom_data_pt(0));
    }
 
    GeomObject** all_geom_object_pt()
    {
      return &(Spine_pt->geom_object_pt(0));
    }
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  class SpineFiniteElement
  {
  public:
    SpineFiniteElement() {}
 
    virtual ~SpineFiniteElement() {}
  };
 
 
  //========================================================================
  //========================================================================
  template<class ELEMENT>
  class SpineElement
    : public ElementWithSpecificMovingNodes<ELEMENT, SpineNode>,
      public SpineFiniteElement
  {
  private:
    unsigned* Spine_geometric_index;
 
    void complete_setup_of_dependencies();
 
 
  public:
    SpineElement()
      : ElementWithSpecificMovingNodes<ELEMENT, SpineNode>(),
        SpineFiniteElement(),
        Spine_geometric_index(0)
    {
    }
 
    SpineElement(FiniteElement* const& element_pt, const int& face_index)
      : ElementWithSpecificMovingNodes<ELEMENT, SpineNode>(element_pt,
                                                           face_index),
        SpineFiniteElement(),
        Spine_geometric_index(0)
    {
    }
 
    ~SpineElement()
    {
      if (Spine_geometric_index)
      {
        delete[] Spine_geometric_index;
      }
    }
 
    inline int spine_local_eqn(const unsigned& n)
    {
#ifdef RANGE_CHECKING
      const unsigned n_node = this->nnode();
      if (n >= n_node)
      {
        std::ostringstream error_message;
        error_message << "Range Error:  Node number " << n
                      << " is not in the range (0," << n_node - 1 << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
#ifdef PARANOID
      // If there is no spine then you can't get the local equation
      if (Spine_geometric_index[n] == this->ngeom_data())
      {
        std::ostringstream error_stream;
        error_stream << "SpineNode " << n
                     << " does not have a Spine attached,\n"
                     << "so you can't get its local equation number.\n"
                     << "Check that the Mesh is correctly associating Spines "
                        "with is Nodes\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return this->geometric_data_local_eqn(Spine_geometric_index[n], 0);
    }
  };
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<SpineElement<ELEMENT>>
    : public virtual FaceGeometry<ELEMENT>
  {
  public:
    FaceGeometry() : FaceGeometry<ELEMENT>() {}
 
  protected:
  };
 
  //=====================================================================
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<FaceGeometry<SpineElement<ELEMENT>>>
    : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT>>() {}
 
  protected:
  };
 
  //=====================================================================
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<SpineElement<FaceGeometry<ELEMENT>>>
    : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT>>() {}
 
  protected:
  };
 
 
 
 
  //========================================================================
  //========================================================================
  class SpineMesh : public virtual Mesh
  {
  protected:
    Vector<Spine*> Spine_pt;
 
  public:
    virtual ~SpineMesh();
 
    Spine*& spine_pt(const unsigned long& i)
    {
      return Spine_pt[i];
    }
 
    const Spine* spine_pt(const unsigned long& i) const
    {
      return Spine_pt[i];
    }
 
    unsigned long nspine() const
    {
      return Spine_pt.size();
    }
 
    void add_spine_pt(Spine* const& spine_pt)
    {
      Spine_pt.push_back(spine_pt);
    }
 
    // Can safely cast the nodes to SpineNodes
    SpineNode* node_pt(const unsigned long& n)
    {
#ifdef PARANOID
      if (!dynamic_cast<SpineNode*>(Node_pt[n]))
      {
        std::ostringstream error_message;
        error_message << "Node " << n << "is a " << typeid(Node_pt[n]).name()
                      << ", not a SpineNode" << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Return a cast to the pointer to the node
      return (dynamic_cast<SpineNode*>(Node_pt[n]));
    }
 
    SpineNode* element_node_pt(const unsigned long& e, const unsigned& n)
    {
#ifdef PARANOID
      // Try to cast to FiniteElement
      FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
      if (el_pt == 0)
      {
        throw OomphLibError(
          "Can't execute element_node_pt(...) for non FiniteElements",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
      if (!dynamic_cast<SpineNode*>(el_pt->node_pt(n)))
      {
        std::ostringstream error_message;
        error_message << "Node " << n << "is a "
                      << typeid(el_pt->node_pt(n)).name() << ", not a SpineNode"
                      << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Return a cast to the node pointer
      return (dynamic_cast<SpineNode*>(
        dynamic_cast<FiniteElement*>(Element_pt[e])->node_pt(n)));
    }
 
    // N.B.: Since SpineElement<ELEMENT>'s are templated, we need the template
    // in this function so that we can do the dynamic cast to a
    // SpineElement<ELEMENT> template<class ELEMENT> void
    // add_spine_to_element(const unsigned long &e, Spine* spine)
    // {dynamic_cast<ELEMENT*>(Element_pt[e])->add_spine(spine);}
 
    unsigned long assign_global_spine_eqn_numbers(Vector<double*>& Dof_pt);
 
    void describe_spine_dofs(std::ostream& out,
                             const std::string& current_string) const;
 
    void set_mesh_level_time_stepper(TimeStepper* const& time_stepper_pt,
                                     const bool& preserve_existing_data)
    {
      this->set_spine_time_stepper(time_stepper_pt, preserve_existing_data);
    }
 
    void set_spine_time_stepper(TimeStepper* const& time_stepper_pt,
                                const bool& preserve_existing_data);
 
    void set_consistent_pinned_spine_values_for_continuation(
      ContinuationStorageScheme* const& continuation_stepper_pt);
 
 
    bool does_pointer_correspond_to_spine_data(double* const& parameter_pt);
 
    void node_update(const bool& update_all_solid_nodes = false);
 
    virtual void spine_node_update(SpineNode* spine_node_pt) = 0;
 
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Woverloaded-virtual"
#endif
 
    void dump(std::ofstream& dump_file) const;
 
#ifdef __clang__
#pragma clang diagnostic pop
#endif
 
    void read(std::ifstream& restart_file);
  };
 
 
  // Functions for the SpineElement class
 
 
  //=================================================================
  //=================================================================
  template<class ELEMENT>
  void SpineElement<ELEMENT>::complete_setup_of_dependencies()
  {
    // Call function of underlying element
    ElementWithSpecificMovingNodes<ELEMENT,
                                   SpineNode>::complete_setup_of_dependencies();
 
    // Sort out the spine index stuff
    // Find the data that correspond to spine heights
    {
      const unsigned n_node = this->nnode();
      // Allocate memory
      if (Spine_geometric_index)
      {
        delete[] Spine_geometric_index;
      }
      Spine_geometric_index = new unsigned[n_node];
 
      // Now loop over data and find out where it fits
      for (unsigned n = 0; n < n_node; n++)
      {
        // Find pointer to the spine
        Spine* const spine_pt =
          static_cast<SpineNode*>(this->node_pt(n))->spine_pt();
 
        // If there is a spine then find the pointer to the data
        if (spine_pt)
        {
          // Find the pointer to the data
          Data* spine_height_data_pt = spine_pt->spine_height_pt();
 
          // Now find the index of the corresponding spine
          const unsigned n_node_update_data = this->ngeom_data();
          for (unsigned i = 0; i < n_node_update_data; i++)
          {
            if (this->Geom_data_pt[i] == spine_height_data_pt)
            {
              Spine_geometric_index[n] = i;
              break;
            }
          }
        }
        // Otherwise issue a warning
        else
        {
          // Set the spine_geometric_index out of range,
          // which will cause the spine_local_eqn to return a pinned value
          Spine_geometric_index[n] = this->ngeom_data();
        }
      }
    }
  }
 
} // namespace oomph
 
#endif