hijacked_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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 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 hijacked elements

//Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_HIJACKED_ELEMENTS_HEADER
#define OOMPH_HIJACKED_ELEMENTS_HEADER


//oomph-lib header
#include "elements.h"
#include "spines.h"

namespace oomph
{

//========================================================================
/// HijackedElement base class that provides storage and access funcitons 
/// for pointers to the global equation numbers that are hijacked by 
/// the HijackedElement. A default residuals multiplier is also provided.
//========================================================================
class HijackedElementBase
{
  public:
 
 /// \short Constructor, initialise the pointer to the equation numbers
 /// for the storage to zero
 HijackedElementBase() : Hijacked_global_eqn_number_pt(0), 
  Hijacked_local_eqn_number_pt(0) 
  {
   //Set the default value of the pointer to the residuals multiplier.
   //The default value is zero, so that the element is "completely hijacked"
   //and the exisiting contributions to the residuals and jacobian are
   //totally wiped
   Residual_multiplier_pt = &Default_residual_multiplier;
  }
 
 /// Destructor, destroy the storage for the equation numbers
 virtual ~HijackedElementBase();

 /// \short Reset the hijacked data pt, so that none of the equations in the 
 /// element are hijacked.
 void unhijack_all_data() {Hijacked_global_eqn_number_pt->clear();}
 
 /// Return the value of the residual multiplier
 inline const double &residual_multiplier() const 
  {return *Residual_multiplier_pt;}
 
 /// Return the pointer to the residual multiplier
 inline double* &residual_multiplier_pt() {return Residual_multiplier_pt;}

  protected:

 /// Pointer to a Set of pointers to the equation numbers that will be
 /// hijacked by this element. Note that these MUST be pointers because
 /// hijacking information is set BEFORE global equation numbers have 
 /// been assigned.
 std::set<long*> *Hijacked_global_eqn_number_pt;

 /// \short Pointer to  a vector of integers 
 /// containing the local equation numbers of 
 /// any hijacked variables in the element.
 Vector<int> *Hijacked_local_eqn_number_pt;
  
 /// \short Pointer to a double that multiplies the contribution to
 /// the residuals from the original element.
 /// This is usually used as a homotopy parameter to permit a smooth
 /// transition between different types of boundary conditions, rather than
 /// switching them on or off abruptly
 double *Residual_multiplier_pt;
  
 /// Static default value for the double that multiplies the original residuals
 static double Default_residual_multiplier;

 /// \short Mark the global equation, addressed by global_eqn_pt, 
 /// as hijacked by this element.
 void hijack_global_eqn(long* const &global_eqn_pt);
 
 /// \short The global equation, addressed by global_eqn_pt,
 /// is no longer hijacked by this element.
 void unhijack_global_eqn(long* const &global_eqn_pt);
 
};


//========================================================================
/// \short Hijacked elements are elements in which one or more
/// Data values that affect the element's residuals, are determined
/// by another element -- the data values are then said to have
/// been hijacked by another element. The main functionality
/// added by the Hijacked element class is that it wipes out
/// those entries in the element's residual vector and those rows in
/// the element's Jacobian matrix that are determined by the "other" 
/// elements that have hijacked the values. Note that for continuation
/// in homotopy parameters, it may be desriable to multiply the residuals 
/// and corresponding jacobian entries by a "homotopy parameter". The
/// value of this parameter can be set by assigning residual_multiplier_pt()
/// which has a default value of zero. Note: it would be possible to extend 
/// the functionality so that different residuals are multiplied by 
/// different values, but will this ever be required?
//========================================================================
template<class ELEMENT>
class Hijacked : public virtual ELEMENT, public virtual HijackedElementBase
{

public:

 /// Constructor, call the constructors of the base elements
 Hijacked() : ELEMENT(), HijackedElementBase() {}

 /// Constructor used for hijacking face elements
 Hijacked(FiniteElement* const &element_pt, const int &face_index) : 
  ELEMENT(element_pt, face_index),
  HijackedElementBase() {}


 /// \short Constructor used for hijacking face elements with specification of
 /// ID of additional variables
  Hijacked(FiniteElement* const &element_pt, 
           const int &face_index, 
           const unsigned &id=0) : 
  ELEMENT(element_pt, face_index,id),
  HijackedElementBase() {}
 
  
 /// \short Hijack the i-th value stored at internal data n.
 /// Optionally return a custom-made (copied) data object that contains only 
 /// the hijacked value. This can be used as the input to other elements.
 /// Note that the calling program assumes responsibility for this 
 /// data object and must clean it up. The default is that
 /// the data object is returned
  Data *hijack_internal_value(const unsigned &n, const unsigned &i,
                              const bool &return_data=true)
  {
   //Initialise pointer to zero
   Data *temp_data_pt = 0;

   //If desired,
   //Create a new Data object containing only the value that is to be hijacked
   if(return_data)
     {temp_data_pt = new HijackedData(i,this->internal_data_pt(n));}
 
   //Mark the value as hijacked
   hijack_global_eqn(this->internal_data_pt(n)->eqn_number_pt(i));

   //Return the pointer to the data
   return temp_data_pt;
  }
 

 /// \short Hijack the i-th value stored at external data n.
 /// Optionally return a custom-made (copied) data object that contains only 
 /// the hijacked value. Note that the calling program assumes
 /// responsibility for this data object and must clean it up.
 /// The default is that the data object is returned
  Data *hijack_external_value(const unsigned &n, const unsigned &i,
                              const bool &return_data=true)
  {
    //Initialise pointer to zero
   Data *temp_data_pt = 0;
   //If desired
   //create a new Data object containing only the value that is to be hijacked
   if(return_data)
     {temp_data_pt = new HijackedData(i,this->external_data_pt(n));}
   
   //Mark the value as hijacked
   hijack_global_eqn(this->external_data_pt(n)->eqn_number_pt(i));
   
   //Return the pointer to the data
   return temp_data_pt;
  }

 /// \short Hijack the i-th value stored at node n.
 /// Optionally return a custom-made (copied) data object that contains only 
 /// the hijacked value. Once again, the calling program must
 /// clean up the allocated Data object.
 /// The default is that the data object is returned
  Data* hijack_nodal_value(const unsigned &n, const unsigned &i,
                           const bool &return_data=true)
  {
    //Initialise pointer to zero
   Data *temp_data_pt = 0;
   //If desired
   //create a new Data object containing only the value that is to be hijacked
   if(return_data)
     {temp_data_pt = new HijackedData(i,this->node_pt(n));}
   
   //Mark the value as hijacked
   hijack_global_eqn(this->node_pt(n)->eqn_number_pt(i));
   
   //Return the pointer to the data
   return temp_data_pt;
  }
 
 /// \short Hijack the i-th positional value stored at node n.
 /// Optionaly return a custom-made (copied) data object that contains only 
 /// the hijacked value. Again, responsibility for the memory allocated
 /// lies with the calling function.
 /// The default is that the data object is returned
 Data* hijack_nodal_position_value(const unsigned &n, 
                                   const unsigned &i,
                                   const bool &return_data=true)
  {
   //Can we do the casting?
   SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(this->node_pt(n));
   if(solid_node_pt == 0)
    {
     std::string error_message = "Failed to cast to SolidNode\n ";
     error_message +=
      "You may be trying to hijack a non-elastic element\n";

     throw OomphLibError(error_message,
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }

    //Initialise pointer to zero
   Data *temp_data_pt = 0;
   //If desired
   //create a new Data object containing only the value that is to be hijacked
   if(return_data)
    {temp_data_pt = 
        new HijackedData(i,solid_node_pt->variable_position_pt());}
   
   //Mark the value as hijacked
   hijack_global_eqn(solid_node_pt->variable_position_pt()->eqn_number_pt(i));

   //Return the pointer to the data
   return temp_data_pt;
  }

 /// \short Hijack the i-th value stored at the spine that affects
 /// local node n.
 /// Optionally return a custom-made (copied) data object that contains only 
 /// the hijacked value. Deletion must be handled at the higher level
 /// The default is that the data object is returned
  Data* hijack_nodal_spine_value(const unsigned &n, const unsigned &i,
                                 const bool &return_data=true)
  {
   //Can we actually do this casting
   SpineNode* spine_node_pt = dynamic_cast<SpineNode*>(this->node_pt(n));
   if(spine_node_pt == 0)
    {
     std::string error_message = "Failed to cast to SpineNode\n ";
     error_message +=
      "You may be trying to hijack a non-spine element\n";

     throw OomphLibError(error_message,
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }

    //Initialise pointer to zero
   Data *temp_data_pt = 0;
   //If desired
   //create a new Data object containing only the value that is to be hijacked
   if(return_data)
    {
     temp_data_pt = new HijackedData(i,spine_node_pt->
                                     spine_pt()->spine_height_pt());
    }
      
   //Mark the data as hijacked
   hijack_global_eqn(spine_node_pt->spine_pt()->
                     spine_height_pt()->eqn_number_pt(i));

   //Return the pointer to the data
   return temp_data_pt;
  }


 /// \short Set up the local equation numbers for the underlying element,
 /// then set up the local arrays to hold the hijacked variables.
 /// If the boolean argument is true then pointers to the associated degrees
 /// of freedom are stored in the array Dof_pt
 void assign_local_eqn_numbers(const bool &store_local_dof_pt)
  {
   //If things have already been allocated,
   //clear the local hijacked array, so that if the equation numbers
   //are reassigned after changes in the boundary conditions, the
   //correct terms will be in the array.
   if(Hijacked_local_eqn_number_pt!=0)
    {Hijacked_local_eqn_number_pt->clear();}
   //Otherwise allocate it
   else
    {Hijacked_local_eqn_number_pt = new Vector<int>;}
    

   //Call the hijacked element's assign_local_eqn_numbers
   ELEMENT::assign_local_eqn_numbers(store_local_dof_pt);
   
   //If any values have been hijacked, we need to find the corresponding
   //local equation numbers
   if(Hijacked_global_eqn_number_pt!=0)
    {
     //Now loop over the local array that stores GLOBAL equation numbers
     //to check if any of the local values have in fact been hijacked
     // by "somebody"
     unsigned n_dof = this->ndof();
     for(unsigned i=0;i<n_dof;i++)
      {
       //Loop over *all* hijacked data
       for(std::set<long*>::iterator it=Hijacked_global_eqn_number_pt->begin();
           it!=Hijacked_global_eqn_number_pt->end();++it)
        {

         //If the hijacked (global!) equation is not pinned
         if(*(*it) >= 0)
          {
           //Get the (unsigned) global equation number: 
           //Note that it is only unsigned AFTER the test above.
           unsigned long hijacked_eqn_number= *(*it);
           
           //If a GLOBAL variable used by this element is hijacked add the 
           //variable's LOCAL index to the array
           if(hijacked_eqn_number == this->eqn_number(i))
          { 
           Hijacked_local_eqn_number_pt->push_back(i);
           break;
          }
          }
        }
      }
    }
  }

 /// \short Get the residuals from the underlying element, but then wipe the
 /// entries in the residual vector that correspond to hijacked
 /// values -- they will be computed by other elements.
 void get_residuals(Vector<double> &residuals)
  {
   //Get parent redisuals
   ELEMENT::get_residuals(residuals);
   //Multiply any hijacked dofs by the residual multiplier
   //(default value is zero)
   unsigned n_hijacked = Hijacked_local_eqn_number_pt->size();
   for(unsigned i=0;i<n_hijacked;i++)
    {
     residuals[(*Hijacked_local_eqn_number_pt)[i]] *= residual_multiplier();
    }
  }


 /// \short Get the residuals and Jacobian matrix from the underlying 
 /// element, but then wipe the entries in the residual vector and the 
 /// rows in the Jacobian matrix that correspond to hijacked
 /// values -- they will be computed by other elements.
 void get_jacobian(Vector<double> &residuals, 
                   DenseMatrix<double> &jacobian)
  {
   //Call the element's get jacobian function
   ELEMENT::get_jacobian(residuals,jacobian);
   //Wipe any hijacked dofs
   unsigned n_hijacked = Hijacked_local_eqn_number_pt->size();
   unsigned n_dof = this->ndof();
   for(unsigned i=0;i<n_hijacked;i++)
    {
     //Multiply any hijacked dofs by the residual multiplier
     //(default value is zero)
     residuals[(*Hijacked_local_eqn_number_pt)[i]] *= residual_multiplier();
     //Multiply the row in the Jacobian matrix by the residual
     //multiplier for consistency
     for(unsigned j=0;j<n_dof;j++)
      {
       jacobian((*Hijacked_local_eqn_number_pt)[i],j) *= residual_multiplier();
      }
    }
  }


};



//============================================================================
///\short Explicit definition of the face geometry of hijacked elements:
///the same as the face geometry of the underlying element
//============================================================================
template<class ELEMENT>
class FaceGeometry<Hijacked<ELEMENT> >:
public virtual FaceGeometry<ELEMENT>  
{
  public:

 /// Constructor
 FaceGeometry() : FaceGeometry<ELEMENT>() {}

  protected:
};

//============================================================================
///\short Explicit definition of the face geometry of hijacked elements:
///the same as the face geometry of the underlying element
//============================================================================
template<class ELEMENT>
class FaceGeometry<FaceGeometry<Hijacked<ELEMENT> > >:
public virtual FaceGeometry<FaceGeometry<ELEMENT> >  
{
  public:

 /// Constructor
 FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT> >() {}

  protected:
};

//============================================================================
///\short Explicit definition of the face geometry of hijacked elements:
///the same as the face geometry of the underlying element
//============================================================================
template<class ELEMENT>
class FaceGeometry<Hijacked<FaceGeometry<ELEMENT> > >:
public virtual FaceGeometry<FaceGeometry<ELEMENT> >  
{
  public:

 /// Constructor
 FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT> >() {}

  protected:
};


}

#endif