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// 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 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
{
  //========================================================================
  //========================================================================
  class HijackedElementBase
  {
  public:
    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;
    }
 
    virtual ~HijackedElementBase();
 
    void unhijack_all_data()
    {
      Hijacked_global_eqn_number_pt->clear();
    }
 
    inline const double& residual_multiplier() const
    {
      return *Residual_multiplier_pt;
    }
 
    inline double*& residual_multiplier_pt()
    {
      return Residual_multiplier_pt;
    }
 
  protected:
    std::set<long*>* Hijacked_global_eqn_number_pt;
 
    Vector<int>* Hijacked_local_eqn_number_pt;
 
    double* Residual_multiplier_pt;
 
    static double Default_residual_multiplier;
 
    void hijack_global_eqn(long* const& global_eqn_pt);
 
    void unhijack_global_eqn(long* const& global_eqn_pt);
  };
 
 
  //========================================================================
  //========================================================================
  template<class ELEMENT>
  class Hijacked : public virtual ELEMENT, public virtual HijackedElementBase
  {
  public:
    Hijacked() : ELEMENT(), HijackedElementBase() {}
 
    Hijacked(FiniteElement* const& element_pt, const int& face_index)
      : ELEMENT(element_pt, face_index), HijackedElementBase()
    {
    }
 
 
    Hijacked(FiniteElement* const& element_pt,
             const int& face_index,
             const unsigned& id = 0)
      : ELEMENT(element_pt, face_index, id), HijackedElementBase()
    {
    }
 
 
    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;
    }
 
 
    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;
    }
 
    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, which may be null
      return temp_data_pt;
    }
 
    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;
    }
 
    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;
    }
 
 
    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;
              }
            }
          }
        }
      }
    }
 
    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();
      }
    }
 
 
    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();
        }
      }
    }
  };
 
 
  //============================================================================
  //============================================================================
  template<class ELEMENT>
  class FaceGeometry<Hijacked<ELEMENT>> : public virtual FaceGeometry<ELEMENT>
  {
  public:
    FaceGeometry() : FaceGeometry<ELEMENT>() {}
 
  protected:
  };
 
  //============================================================================
  //============================================================================
  template<class ELEMENT>
  class FaceGeometry<FaceGeometry<Hijacked<ELEMENT>>>
    : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT>>() {}
 
  protected:
  };
 
  //============================================================================
  //============================================================================
  template<class ELEMENT>
  class FaceGeometry<Hijacked<FaceGeometry<ELEMENT>>>
    : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
  {
  public:
    FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT>>() {}
 
  protected:
  };
 
 
} // namespace oomph
 
#endif