perturbed_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 perturbed spine nodes, elements and meshes
 
#ifndef OOMPH_PERTURBED_SPINES_HEADER
#define OOMPH_PERTURBED_SPINES_HEADER
 
// oomph-lib includes
#include "generic.h"
 
namespace oomph
{
 
 // Forward declaration
 class PerturbedSpineNode;
 
//=======================================================================
//=======================================================================
class PerturbedSpine
{
  
  public:
  
 PerturbedSpine(Spine*& base_spine_pt) : Base_spine_pt(base_spine_pt)
  {
   Geom_data_pt.resize(1);
   
   // Create Data for the two "heights". By default they are free
   Geom_data_pt[0] = new Data(2);
  }
  
 PerturbedSpine(TimeStepper* const &time_stepper_pt,
                Spine*& base_spine_pt) : Base_spine_pt(base_spine_pt)
  {
   Geom_data_pt.resize(1);
   
   // Create Data for the two "heights". By default they are free
   Geom_data_pt[0] = new Data(time_stepper_pt,2);
  }
    
 PerturbedSpine(Spine*& base_spine_pt,
                const double& height) : Base_spine_pt(base_spine_pt)
  {
   Geom_data_pt.resize(1);
   
   // Create Data for the two "heights". By default they are free
   Geom_data_pt[0] = new Data(2);
   
   // Set value
   for(unsigned i=0;i<2;i++) { Geom_data_pt[0]->set_value(i,height); }
  }
     
 ~PerturbedSpine() { delete Geom_data_pt[0]; }
 
 Spine* &base_spine_pt() { return Base_spine_pt; }
 
 double& height(const unsigned& i)
  {
   return *(Geom_data_pt[0]->value_pt(i));
  }
 
 double& height(const unsigned& t, const unsigned& i)
  {
   return *(Geom_data_pt[0]->value_pt(t,i));
  }
 
 Data*& height_pt() { return Geom_data_pt[0]; }
 
 Data* height_pt() const { return Geom_data_pt[0]; }
 
 PerturbedSpineNode*& node_at_bottom_of_spine_pt()
  {
   return Node_at_bottom_of_spine_pt;
  }
 
 PerturbedSpineNode* node_at_bottom_of_spine_pt() const
  {
   return Node_at_bottom_of_spine_pt;
  }
 
 PerturbedSpineNode*& node_at_top_of_spine_pt()
  {
   return Node_at_top_of_spine_pt;
  }
 
 PerturbedSpineNode* node_at_top_of_spine_pt() const
  {
   return Node_at_top_of_spine_pt;
  }
 
 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:
 
 Spine* Base_spine_pt;
 
 Vector<Data*> Geom_data_pt;
 
 Vector<GeomObject*> Geom_object_pt;
 
 Vector<double> Geom_parameter;
 
 PerturbedSpineNode* Node_at_bottom_of_spine_pt;
 PerturbedSpineNode* Node_at_top_of_spine_pt;
 
}; // End of PerturbedSpine class definition
 
 
 
// Forward declaration
class PerturbedSpineMesh;
 
//=======================================================================
//=======================================================================
class PerturbedSpineNode : public Node
{
  private:
 
 PerturbedSpine* PerturbedSpine_pt;
 
 double Fraction;
 
 PerturbedSpineMesh* PerturbedSpine_mesh_pt;
 
 unsigned Node_update_fct_id;
 
 
  public:
 
 PerturbedSpineNode(const unsigned &n_dim,
                    const unsigned &n_position_type,
                    const unsigned &initial_nvalue) :
  Node(n_dim,n_position_type,initial_nvalue),
  PerturbedSpine_pt(0), Fraction(0),
  PerturbedSpine_mesh_pt(0), Node_update_fct_id(0) {}
  
 PerturbedSpineNode(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),
  PerturbedSpine_pt(0),
  Fraction(0), PerturbedSpine_mesh_pt(0), Node_update_fct_id(0) {}
   
 PerturbedSpine* &perturbed_spine_pt() { return PerturbedSpine_pt; }
   
 double &fraction() { return Fraction; }
   
 unsigned& node_update_fct_id() { return Node_update_fct_id; }
   
 PerturbedSpineMesh*& spine_mesh_pt() { return PerturbedSpine_mesh_pt; }
 
 double &h(const unsigned& i) { return PerturbedSpine_pt->height(i); }
   
 void node_update(const bool& update_all_time_levels_for_new_node=false);
   
 unsigned ngeom_data() const
 {
  if(PerturbedSpine_pt) { return PerturbedSpine_pt->ngeom_data(); }
  else { return 0; }
 }
 
 unsigned ngeom_object() const
 {
  if(PerturbedSpine_pt) { return PerturbedSpine_pt->ngeom_object(); }
  else { return 0; }
 }
 
 Data** all_geom_data_pt() { return &(PerturbedSpine_pt->geom_data_pt(0)); }
   
 GeomObject** all_geom_object_pt()
  {
   return &(PerturbedSpine_pt->geom_object_pt(0));
  }
   
}; // End of PerturbedSpineNode class definition
 
 
 
 
 
 
//=======================================================================
//=======================================================================
class PerturbedSpineFiniteElement
{
  public:
 
 PerturbedSpineFiniteElement() {}
 
 virtual ~PerturbedSpineFiniteElement(){}
 
};
 
 
 
//=======================================================================
//=======================================================================
template<class ELEMENT>
class PerturbedSpineElement : 
public ElementWithSpecificMovingNodes<ELEMENT,PerturbedSpineNode>,
 public PerturbedSpineFiniteElement
{
  private:
 
 unsigned *PerturbedSpine_geometric_index;
 
 void complete_setup_of_dependencies();
 
 
  public:
 
 PerturbedSpineElement() : 
  ElementWithSpecificMovingNodes<ELEMENT,PerturbedSpineNode>(),  
  PerturbedSpineFiniteElement(),
  PerturbedSpine_geometric_index(0) {}
 
 PerturbedSpineElement(FiniteElement* const &element_pt,
                       const int &face_index) : 
  ElementWithSpecificMovingNodes<ELEMENT,PerturbedSpineNode>(element_pt,
                                                             face_index),
  PerturbedSpineFiniteElement(),
  PerturbedSpine_geometric_index(0) {}
 
 ~PerturbedSpineElement()
  {
   if(PerturbedSpine_geometric_index) 
    {
     delete[] PerturbedSpine_geometric_index;
    }
  }
 
 inline int spine_local_eqn(const unsigned &n, const unsigned &i)
  {
#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(),
                         "PerturbedSpineElement::spine_local_eqn()",
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
#ifdef PARANOID
   // The perturbed spine "height" only has two components
   if(i>=2)
    {
     std::ostringstream error_message;
     error_message << "Perturbed spines only have two height components:"
                   << "\nComponent i=" << i
                   << " is not in the range (0,2)";
     throw OomphLibError(error_message.str(),
                         "PerturbedSpineElement::spine_local_eqn()",
                         OOMPH_EXCEPTION_LOCATION);
    }
   // If there is no spine then you can't get the local equation
   if(PerturbedSpine_geometric_index[n]==this->ngeom_data())
    {
     std::ostringstream error_stream;
     error_stream 
      << "PerturbedSpineNode " << n << " does not have a PerturbedSpine\n"
      << "attached, so you can't get its local equation number. Check that\n"
      << "the Mesh is correctly associating PerturbedSpines with is Nodes\n";
     throw OomphLibError(error_stream.str(),
                         "PerturbedSpineElement<ELEMENT>::spine_local_eqn()",
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   return this->geometric_data_local_eqn(PerturbedSpine_geometric_index[n],i);
  }
 
 
 virtual int get_local_eqn_number_corresponding_to_geometric_dofs(
  const unsigned &n, const unsigned &i)
 {
  return this->spine_local_eqn(n,i);
 }
 
}; // End of PerturbedSpineElement<ELEMENT> class definition
 
 
 
//=======================================================================
//=======================================================================
template<class ELEMENT>
class FaceGeometry<PerturbedSpineElement<ELEMENT> >:
public virtual FaceGeometry<ELEMENT>  
{
  public:
 
 FaceGeometry() : FaceGeometry<ELEMENT>() {}
};
 
//=======================================================================
//=======================================================================
template<class ELEMENT>
class FaceGeometry<FaceGeometry<PerturbedSpineElement<ELEMENT> > >:
public virtual FaceGeometry<FaceGeometry<ELEMENT> > 
{
  public:
 
 FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT> >() {}
};
 
//=======================================================================
//=======================================================================
template<class ELEMENT>
class FaceGeometry<PerturbedSpineElement<FaceGeometry<ELEMENT> > >:
public virtual FaceGeometry<FaceGeometry<ELEMENT> > 
{
  public:
 
 FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT> >() {}
};
 
 
 
 
 
 
//========================================================================
//========================================================================
class PerturbedSpineMesh : public virtual Mesh
{
  protected:
 
 Vector<PerturbedSpine*> PerturbedSpine_pt;
 
  public:
 
 virtual ~PerturbedSpineMesh()
  {
   // Set the range of PerturbedSpine_pt
   unsigned long PerturbedSpine_pt_range = PerturbedSpine_pt.size();
 
   // Loop over the entries in reverse and free memory
   for(unsigned long i=PerturbedSpine_pt_range;i>0;i--)
    {
     delete PerturbedSpine_pt[i-1]; PerturbedSpine_pt[i-1] = 0;
    }
  }
 
 PerturbedSpine* const & perturbed_spine_pt(const unsigned long &i) const
  {
   return PerturbedSpine_pt[i];
  }
 
 unsigned long nspine() const { return PerturbedSpine_pt.size(); }
 
 void add_perturbed_spine_pt(PerturbedSpine* const &spine_pt)
 {
  PerturbedSpine_pt.push_back(spine_pt);
 }
 
 PerturbedSpineNode* node_pt(const unsigned long &n) 
  {
#ifdef PARANOID
   if(!dynamic_cast<PerturbedSpineNode*>(Node_pt[n]))
    {
     std::ostringstream error_message;
     error_message << "Node " << n << "is a "
                   << typeid(Node_pt[n]).name() 
                   << ", not a PerturbedSpineNode" << std::endl;
     
     throw OomphLibError(error_message.str(),
                         "PerturbedSpineMesh::node_pt()",
                         OOMPH_EXCEPTION_LOCATION);
    } 
#endif
 
   // Return a cast to the pointer to the node
   return (dynamic_cast<PerturbedSpineNode*>(Node_pt[n]));
  }
 
 PerturbedSpineNode* 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",
      "PerturbedSpineMesh::element_node_pt()",
      OOMPH_EXCEPTION_LOCATION);
    }
   if(!dynamic_cast<PerturbedSpineNode*>(el_pt->node_pt(n)))
    {
     std::ostringstream error_message;
     error_message << "Node " << n << "is a "
                   << typeid(el_pt->node_pt(n)).name() 
                   << ", not a PerturbedSpineNode" << std::endl;
     
     throw OomphLibError(error_message.str(),
                         "PerturbedSpineMesh::node_pt()",
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   // Return a cast to the node pointer
   return(dynamic_cast<PerturbedSpineNode*>(
           dynamic_cast<FiniteElement*>(Element_pt[e])->node_pt(n)));
  }
 
 unsigned long assign_global_eqn_numbers(Vector<double *> &Dof_pt);
 
 void node_update(const bool& update_all_solid_nodes=false);
 
 virtual void perturbed_spine_node_update(PerturbedSpineNode* spine_node_pt)=0;
 
 void dump(std::ofstream &dump_file,const bool &use_old_ordering=true) const;
 
 void read(std::ifstream &restart_file);
 
}; // End of PerturbedSpineMesh class definition
 
 
 
 
} // End of oomph namespace
 
#endif