element_with_moving_nodes.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 base class that specified common interfaces
// used in elements with moving Nodes (Spine,Algebraic,MacroElementNodeUpdate)
 
// Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_ELEMENT_WITH_MOVING_NODES
#define OOMPH_ELEMENT_WITH_MOVING_NODES
 
 
// oomph-lib headers
#include "elements.h"
 
namespace oomph
{
 
 
  //=======================================================================
  //=======================================================================
  class ElementWithMovingNodes : public virtual FiniteElement
  {
  public:
    enum
    {
      Shape_derivs_by_chain_rule,
      Shape_derivs_by_direct_fd,
      Shape_derivs_by_fastest_method
    };
 
    ElementWithMovingNodes()
      : Geometric_data_local_eqn(0),
        Bypass_fill_in_jacobian_from_geometric_data(false),
        Evaluate_dresidual_dnodal_coordinates_by_fd(false),
        Method_for_shape_derivs(Shape_derivs_by_direct_fd)
    // hierher: Anything other than the fd-based method is currently broken;
    // at least for refineable elements -- this all needs to be checked
    // VERY carefully again (see instructions in commit log). Until this
    // has all been fixed/re-validated, we've frozen the default assignment
    // (by breaking the functions that change the setting). Ultimately,
    // we should use:
    // Method_for_shape_derivs(Shape_derivs_by_fastest_method)
    {
    }
 
    ElementWithMovingNodes(const ElementWithMovingNodes&) = delete;
 
    void operator=(const ElementWithMovingNodes&) = delete;
 
    void describe_local_dofs(std::ostream& out,
                             const std::string& current_string) const;
 
    virtual ~ElementWithMovingNodes()
    {
      if (Geometric_data_local_eqn)
      {
        delete[] Geometric_data_local_eqn[0];
        delete[] Geometric_data_local_eqn;
      }
    }
 
    unsigned ngeom_dof() const
    {
      return Ngeom_dof;
    }
 
    inline int geometric_data_local_eqn(const unsigned& n, const unsigned& i)
    {
#ifdef RANGE_CHECKING
      unsigned n_data = Geom_data_pt.size();
      if (n >= n_data)
      {
        std::ostringstream error_message;
        error_message << "Range Error:  Data number " << n
                      << " is not in the range (0," << n_data - 1 << ")";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else
      {
        unsigned n_value = Geom_data_pt[n]->nvalue();
        if (i >= n_value)
        {
          std::ostringstream error_message;
          error_message << "Range Error: value " << i << " at data " << n
                        << " is not in the range (0," << n_value - 1 << ")";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
#endif
#ifdef PARANOID
      // Check that the equations have been allocated
      if (Geometric_data_local_eqn == 0)
      {
        throw OomphLibError(
          "Geometric data local equation numbers have not been allocated",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Geometric_data_local_eqn[n][i];
    }
 
 
    void assemble_set_of_all_geometric_data(
      std::set<Data*>& unique_geom_data_pt);
 
    void identify_geometric_data(std::set<Data*>& geometric_data_pt)
    {
      // Loop over the node update data and add to the set
      const unsigned n_geom_data = Geom_data_pt.size();
      for (unsigned n = 0; n < n_geom_data; n++)
      {
        geometric_data_pt.insert(Geom_data_pt[n]);
      }
    }
 
    bool are_dresidual_dnodal_coordinates_always_evaluated_by_fd() const
    {
      return Evaluate_dresidual_dnodal_coordinates_by_fd;
    }
 
    void enable_always_evaluate_dresidual_dnodal_coordinates_by_fd()
    {
      Evaluate_dresidual_dnodal_coordinates_by_fd = true;
    }
 
    void disable_always_evaluate_dresidual_dnodal_coordinates_by_fd()
    {
      Evaluate_dresidual_dnodal_coordinates_by_fd = false;
    }
 
    void evaluate_shape_derivs_by_direct_fd()
    {
      // hierher: Harmless; simply re-sets the current (and currently only
      // legal) default
      Method_for_shape_derivs = Shape_derivs_by_direct_fd;
    }
 
    void evaluate_shape_derivs_by_chain_rule(
      const bool& i_know_what_i_am_doing = false)
    {
      if (!i_know_what_i_am_doing)
      {
        std::ostringstream error_message;
        error_message
          << "Evaluation of shape derivatives by chain rule is currently \n"
          << "disabled because it's broken, at least for refineable \n"
          << "elements. This all needs to be checked again very carefully\n"
          << "following the instructions in the commit log\n"
          << "If you know what you're doing and want to force this "
             "methodology\n"
          << "call this function with the optional boolean set to true.\n";
        throw OomphLibError(
          error_message.str(),
          "ElementWithMovingNodes::evaluate_shape_derivs_by_chain_rule()",
          OOMPH_EXCEPTION_LOCATION);
      }
      Method_for_shape_derivs = Shape_derivs_by_chain_rule;
    }
 
    void evaluate_shape_derivs_by_fastest_method(
      const bool& i_know_what_i_am_doing = false)
    {
      if (!i_know_what_i_am_doing)
      {
        std::ostringstream error_message;
        error_message
          << "Evaluation of shape derivatives by fastest method is currently \n"
          << "disabled because it's broken, at least for refineable \n"
          << "elements. This all needs to be checked again very carefully\n"
          << "following the instructions in the commit log\n"
          << "If you know what you're doing and want to force this "
             "methodology\n"
          << "call this function with the optional boolean set to true.\n";
        throw OomphLibError(
          error_message.str(),
          "ElementWithMovingNodes::evaluate_shape_derivs_by_fastest_method()",
          OOMPH_EXCEPTION_LOCATION);
      }
      Method_for_shape_derivs = Shape_derivs_by_fastest_method;
    }
 
    int& method_for_shape_derivs()
    {
      return Method_for_shape_derivs;
    }
 
    void enable_bypass_fill_in_jacobian_from_geometric_data()
    {
      Bypass_fill_in_jacobian_from_geometric_data = true;
    }
 
    void disable_bypass_fill_in_jacobian_from_geometric_data()
    {
      Bypass_fill_in_jacobian_from_geometric_data = false;
    }
 
    bool is_fill_in_jacobian_from_geometric_data_bypassed() const
    {
      return Bypass_fill_in_jacobian_from_geometric_data;
    }
 
  protected:
    virtual void get_dnodal_coordinates_dgeom_dofs(
      RankThreeTensor<double>& dnodal_coordinates_dgeom_dofs);
 
    void complete_setup_of_dependencies();
 
    virtual void assign_all_generic_local_eqn_numbers(
      const bool& store_local_dof_pt);
 
    void fill_in_jacobian_from_geometric_data(Vector<double>& residuals,
                                              DenseMatrix<double>& jacobian);
 
    void fill_in_jacobian_from_geometric_data(DenseMatrix<double>& jacobian)
    {
      if (!Bypass_fill_in_jacobian_from_geometric_data)
      {
        // Allocate storage for the residuals
        const unsigned n_dof = ndof();
        Vector<double> residuals(n_dof);
 
        // Get the residuals for the entire element
        get_residuals(residuals);
 
        // Call the jacobian calculation
        fill_in_jacobian_from_geometric_data(residuals, jacobian);
      }
    }
 
 
    Vector<Data*> Geom_data_pt;
 
  public:
    unsigned ngeom_data() const
    {
      return Geom_data_pt.size();
    }
 
  private:
    int** Geometric_data_local_eqn;
 
    unsigned Ngeom_dof;
 
    bool Bypass_fill_in_jacobian_from_geometric_data;
 
    bool Evaluate_dresidual_dnodal_coordinates_by_fd;
 
    int Method_for_shape_derivs;
  };
 
 
  //===============================================================
  //==============================================================
  template<class ELEMENT, class NODE_TYPE>
  class ElementWithSpecificMovingNodes : public ELEMENT,
                                         public ElementWithMovingNodes
  {
  public:
    void describe_local_dofs(std::ostream& out,
                             const std::string& current_string) const
    {
      ELEMENT::describe_local_dofs(out, current_string);
      ElementWithMovingNodes::describe_local_dofs(out, current_string);
    }
 
    ElementWithSpecificMovingNodes() : ELEMENT(), ElementWithMovingNodes() {}
 
    ElementWithSpecificMovingNodes(FiniteElement* const& element_pt,
                                   const int& face_index)
      : ELEMENT(element_pt, face_index), ElementWithMovingNodes()
    {
    }
 
    ~ElementWithSpecificMovingNodes() {}
 
    void describe_local_dofs(std::ostream& out, std::string& curr_str)
    {
      ElementWithMovingNodes::describe_local_dofs(out, curr_str);
      ELEMENT::describe_local_dofs(out, curr_str);
    }
 
 
    Node* construct_node(const unsigned& n)
    {
      // Assign a node to the local node pointer
      // The dimension and number of values are taken from internal element data
      // The number of timesteps to be stored comes from the problem!
      this->node_pt(n) = new NODE_TYPE(this->nodal_dimension(),
                                       this->nnodal_position_type(),
                                       this->required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can find it
      return this->node_pt(n);
    }
 
    Node* construct_node(const unsigned& n, TimeStepper* const& time_stepper_pt)
    {
      // Assign a node to the local node pointer
      // The dimension and number of values are taken from internal element data
      // The number of timesteps to be stored comes from the problem!
      this->node_pt(n) = new NODE_TYPE(time_stepper_pt,
                                       this->nodal_dimension(),
                                       this->nnodal_position_type(),
                                       this->required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can find it
      return this->node_pt(n);
    }
 
    Node* construct_boundary_node(const unsigned& n)
    {
      // Assign a node to the local node pointer
      // The dimension and number of values are taken from internal element data
      // The number of timesteps to be stored comes from the problem!
      this->node_pt(n) =
        new BoundaryNode<NODE_TYPE>(this->nodal_dimension(),
                                    this->nnodal_position_type(),
                                    this->required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can find it
      return this->node_pt(n);
    }
 
    Node* construct_boundary_node(const unsigned& n,
                                  TimeStepper* const& time_stepper_pt)
    {
      // Assign a node to the local node pointer
      // The dimension and number of values are taken from internal element data
      // The number of timesteps to be stored comes from the problem!
      this->node_pt(n) =
        new BoundaryNode<NODE_TYPE>(time_stepper_pt,
                                    this->nodal_dimension(),
                                    this->nnodal_position_type(),
                                    this->required_nvalue(n));
      // Now return a pointer to the node, so that the mesh can find it
      return this->node_pt(n);
    }
 
 
    void complete_setup_of_dependencies()
    {
      // Call function of underlying element
      ELEMENT::complete_setup_of_dependencies();
 
      // Call function of the element with moving nodes
      ElementWithMovingNodes::complete_setup_of_dependencies();
    }
 
 
    void assign_all_generic_local_eqn_numbers(const bool& store_local_dof_pt)
    {
      // Call the generic local equation numbering scheme of the ELEMENT
      ELEMENT::assign_all_generic_local_eqn_numbers(store_local_dof_pt);
      ElementWithMovingNodes::assign_all_generic_local_eqn_numbers(
        store_local_dof_pt);
    }
 
    void get_jacobian(Vector<double>& residuals, DenseMatrix<double>& jacobian)
    {
      ELEMENT::get_jacobian(residuals, jacobian);
 
      // Now call the additional geometric Jacobian terms
      this->fill_in_jacobian_from_geometric_data(jacobian);
    }
 
    void get_jacobian_and_mass_matrix(Vector<double>& residuals,
                                      DenseMatrix<double>& jacobian,
                                      DenseMatrix<double>& mass_matrix)
    {
      ELEMENT::get_jacobian_and_mass_matrix(residuals, jacobian, mass_matrix);
 
      // Now call the additional geometric Jacobian terms
      this->fill_in_jacobian_from_geometric_data(jacobian);
    }
 
 
  private:
  };
 
 
} // namespace oomph
#endif