element_with_external_element.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 the class
// ElementWithExternalElement which stores pointers to external elements
 
// Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_ELEMENT_WITH_EXTERNAL_ELEMENT_HEADER
#define OOMPH_ELEMENT_WITH_EXTERNAL_ELEMENT_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// Oomph-lib headers
#include "elements.h"
 
namespace oomph
{
  //========================================================================
  //========================================================================
  class ElementWithExternalElement : public virtual FiniteElement
  {
  public:
    ElementWithExternalElement()
      : FiniteElement(),
        Add_external_interaction_data(true),
        Add_external_geometric_data(true),
        External_interaction_field_data_pt(0),
        External_interaction_geometric_data_pt(0),
        Ninteraction(0),
        Nintpt(0),
        Nexternal_element_storage(0),
        Nexternal_interaction_field_data(0),
        Nexternal_interaction_geometric_data(0),
        External_element_pt(0),
        External_element_local_coord(0),
        External_interaction_field_data_index(0),
        External_interaction_field_data_local_eqn(0),
        External_interaction_geometric_data_index(0),
        External_interaction_geometric_data_local_eqn(0)
    {
    }
 
    virtual ~ElementWithExternalElement();
 
    ElementWithExternalElement(const ElementWithExternalElement&) = delete;
 
    void operator=(const ElementWithExternalElement&) = delete;
 
    bool storage_has_been_allocated() const
    {
      // If either of the storage arrays is zero, then storage has not
      // yet been allocated. Both arrays are allocated at once, so
      // if one is zero both will (should) be
      if ((External_element_pt == 0) || (External_element_local_coord == 0))
      {
        return false;
      }
      else
      {
        return true;
      }
    }
 
    FiniteElement*& external_element_pt(const unsigned& interaction_index,
                                        const unsigned& ipt)
    {
#ifdef PARANOID
      check_storage_allocated();
#endif
#ifdef RANGE_CHECKING
      range_check(interaction_index, ipt);
#endif
      // Return the appropriate entry in the storage array
      return External_element_pt[Nintpt * interaction_index + ipt];
    }
 
    FiniteElement* const& external_element_pt(const unsigned& interaction_index,
                                              const unsigned& ipt) const
    {
#ifdef PARANOID
      check_storage_allocated();
#endif
#ifdef RANGE_CHECKING
      range_check(interaction_index, ipt);
#endif
      // Return the appropriate entry in the storage array
      return External_element_pt[Nintpt * interaction_index + ipt];
    }
 
    Vector<double>& external_element_local_coord(
      const unsigned& interaction_index, const unsigned& ipt)
    {
#ifdef PARANOID
      check_storage_allocated();
#endif
#ifdef RANGE_CHECKING
      range_check(interaction_index, ipt);
#endif
      return External_element_local_coord[Nintpt * interaction_index + ipt];
    }
 
    Vector<double> const& external_element_local_coord(
      const unsigned& interaction_index, const unsigned& ipt) const
    {
#ifdef PARANOID
      check_storage_allocated();
#endif
#ifdef RANGE_CHECKING
      range_check(interaction_index, ipt);
#endif
      return External_element_local_coord[Nintpt * interaction_index + ipt];
    }
 
 
    void output_external_elements(std::ostream& outfile,
                                  const unsigned& interaction_index);
 
    void initialise_external_element_storage();
 
    void flush_all_external_element_storage();
 
    inline void set_ninteraction(const unsigned& n_interaction)
    {
      Ninteraction = n_interaction;
    }
 
    unsigned ninteraction() const
    {
      return Ninteraction;
    }
 
    virtual void identify_all_field_data_for_external_interaction(
      Vector<std::set<FiniteElement*>> const& external_elements_pt,
      std::set<std::pair<Data*, unsigned>>& paired_interaction_data);
 
    virtual void identify_all_geometric_data_for_external_interaction(
      Vector<std::set<FiniteElement*>> const& external_elements_pt,
      std::set<Data*>& external_geometric_data_pt);
 
    unsigned nexternal_interaction_field_data() const
    {
      return Nexternal_interaction_field_data;
    }
 
    Vector<Data*> external_interaction_field_data_pt() const
    {
      // Find out how many external field data there are
      const unsigned n_external_interaction_field_data =
        nexternal_interaction_field_data();
      // Create and populate  a temporary vector
      Vector<Data*> temp_data(n_external_interaction_field_data);
      for (unsigned i = 0; i < n_external_interaction_field_data; i++)
      {
        temp_data[i] = External_interaction_field_data_pt[i];
      }
      // Return the temporary data
      return temp_data;
    }
 
 
    unsigned nexternal_interaction_geometric_data() const
    {
      return Nexternal_interaction_geometric_data;
    }
 
    bool add_external_geometric_data()
    {
      return Add_external_geometric_data;
    }
 
    bool add_external_interaction_data()
    {
      return Add_external_interaction_data;
    }
 
    Vector<Data*> external_interaction_geometric_data_pt() const
    {
      // Find out how many external field data there are
      const unsigned n_external_interaction_geometric_data =
        nexternal_interaction_geometric_data();
      // Create and populate  a temporary vector
      Vector<Data*> temp_data(n_external_interaction_geometric_data);
      for (unsigned i = 0; i < n_external_interaction_geometric_data; i++)
      {
        temp_data[i] = External_interaction_geometric_data_pt[i];
      }
      // Return the temporary data
      return temp_data;
    }
 
    void ignore_external_geometric_data()
    {
      Add_external_geometric_data = false;
    }
 
    void ignore_external_interaction_data()
    {
      Add_external_interaction_data = false;
    }
 
    void include_external_geometric_data()
    {
      Add_external_geometric_data = true;
    }
 
    void include_external_interaction_data()
    {
      Add_external_interaction_data = true;
    }
 
    bool external_geometric_data_is_included() const
    {
      return Add_external_geometric_data;
    }
 
    void describe_local_dofs(std::ostream& out,
                             const std::string& curr_string) const;
 
  protected:
    void assign_internal_and_external_local_eqn_numbers(
      const bool& store_local_dof_pt)
    {
      // Call the external stuff first so that it is near the front of the
      // list for fast searches when using indexing by global dof for
      // analytic calculation of interaction blocks.
      this->assign_external_interaction_data_local_eqn_numbers(
        store_local_dof_pt);
      // Now call the underlying equation numbering
      GeneralisedElement::assign_internal_and_external_local_eqn_numbers(
        store_local_dof_pt);
    }
 
    void assign_external_interaction_data_local_eqn_numbers(
      const bool& store_local_dof_pt);
 
 
    void fill_in_jacobian_from_external_interaction_by_fd(
      Vector<double>& residuals, DenseMatrix<double>& jacobian)
    {
      // Get the contribution from the external field data
      fill_in_jacobian_from_external_interaction_field_by_fd(residuals,
                                                             jacobian);
      // Get the contribution from the external geometric data
      fill_in_jacobian_from_external_interaction_geometric_by_fd(residuals,
                                                                 jacobian);
    }
 
    void fill_in_jacobian_from_external_interaction_by_fd(
      DenseMatrix<double>& jacobian)
    {
      const unsigned n_dof = this->ndof();
      // Allocate storage for the residuals
      Vector<double> residuals(n_dof, 0.0);
      // Get the residuals for the entire element
      get_residuals(residuals);
      // Call the jacobian calculation
      fill_in_jacobian_from_external_interaction_by_fd(residuals, jacobian);
    }
 
 
    void fill_in_jacobian_from_external_interaction_field_by_fd(
      Vector<double>& residuals, DenseMatrix<double>& jacobian);
 
    void fill_in_jacobian_from_external_interaction_field_by_fd(
      DenseMatrix<double>& jacobian)
    {
      const unsigned n_dof = this->ndof();
      // Allocate storage for the residuals
      Vector<double> residuals(n_dof, 0.0);
      // Get the residuals for the entire element
      get_residuals(residuals);
      // Call the jacobian calculation
      fill_in_jacobian_from_external_interaction_field_by_fd(residuals,
                                                             jacobian);
    }
 
    void fill_in_jacobian_from_external_interaction_geometric_by_fd(
      Vector<double>& residuals, DenseMatrix<double>& jacobian);
 
    void fill_in_jacobian_from_external_interaction_geometric_by_fd(
      DenseMatrix<double>& jacobian)
    {
      const unsigned n_dof = this->ndof();
      // Allocate storage for the residuals
      Vector<double> residuals(n_dof, 0.0);
      // Get the residuals for the entire element
      get_residuals(residuals);
      // Call the jacobian calculation
      fill_in_jacobian_from_external_interaction_geometric_by_fd(residuals,
                                                                 jacobian);
    }
 
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Add the contribution to the residuals
      fill_in_contribution_to_residuals(residuals);
      // Allocate storage for the full residuals  (residuals of entire element)
      const unsigned n_dof = this->ndof();
      Vector<double> full_residuals(n_dof);
      // Get the residuals for the entire element
      get_residuals(full_residuals);
      // Add the internal and external by finite differences
      fill_in_jacobian_from_internal_by_fd(full_residuals, jacobian);
      fill_in_jacobian_from_external_by_fd(full_residuals, jacobian);
      fill_in_jacobian_from_nodal_by_fd(full_residuals, jacobian);
      fill_in_jacobian_from_external_interaction_by_fd(full_residuals,
                                                       jacobian);
    }
 
 
    virtual inline void update_before_external_interaction_field_fd() {}
 
    virtual inline void reset_after_external_interaction_field_fd() {}
 
    virtual inline void update_in_external_interaction_field_fd(
      const unsigned& i)
    {
    }
 
    virtual inline void reset_in_external_interaction_field_fd(
      const unsigned& i)
    {
      update_in_external_interaction_field_fd(i);
    }
 
 
    virtual inline void update_before_external_interaction_geometric_fd() {}
 
    virtual inline void reset_after_external_interaction_geometric_fd() {}
 
    virtual inline void update_in_external_interaction_geometric_fd(
      const unsigned& i)
    {
    }
 
    virtual inline void reset_in_external_interaction_geometric_fd(
      const unsigned& i)
    {
      update_in_external_interaction_geometric_fd(i);
    }
 
    bool Add_external_interaction_data;
 
    bool Add_external_geometric_data;
 
    Data** External_interaction_field_data_pt;
 
    Data** External_interaction_geometric_data_pt;
 
  private:
    void check_storage_allocated() const
    {
      // If either of the storage arrays is zero, then storage has not
      // yet been allocated. Both arrays are allocated at once, so
      // if one is zero both will (should) be
      if ((External_element_pt == 0) || (External_element_local_coord == 0))
      {
        std::ostringstream error_stream;
        error_stream
          << "Storage for the external elements has not been allocated.\n"
          << "initialise_external_element_storage() must be called.\n"
          << "This tends to be done automatically during the multi-domain\n"
          << "setup procedures -- you're most likely to get this error\n"
          << "because you have not specified the number of interactions\n"
          << "that your ElementWithExternalElement is involved in.\n"
          << "You ought to specify this with a call to\n\n"
          << "     ElementWithExternalElement::set_ninteraction(...)\n\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    void range_check(const unsigned& interaction_index,
                     const unsigned& ipt) const
    {
      // Boolean flag used to indicate range error
      bool range_error = false;
      // String for the error message
      std::ostringstream error_message;
      // If there is a range error in the interaction index
      if (interaction_index >= Ninteraction)
      {
        error_message << "Range Error: Interaction " << interaction_index
                      << " is not in the range (0," << Ninteraction - 1 << ")";
        // There has been a range error
        range_error = true;
      }
 
      // If there is a range error in the integration point
      if (ipt >= Nintpt)
      {
        error_message << "Range Error: Integration point " << ipt
                      << " is not in the range (0," << Nintpt - 1 << ")";
        range_error = true;
      }
 
 
      // If there has been a range error report it
      if (range_error)
      {
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    unsigned Ninteraction;
 
    unsigned Nintpt;
 
    unsigned Nexternal_element_storage;
 
    unsigned Nexternal_interaction_field_data;
 
    unsigned Nexternal_interaction_geometric_data;
 
 
    FiniteElement** External_element_pt;
 
    Vector<double>* External_element_local_coord;
 
    unsigned* External_interaction_field_data_index;
 
    int* External_interaction_field_data_local_eqn;
 
    unsigned* External_interaction_geometric_data_index;
 
    int* External_interaction_geometric_data_local_eqn;
  };
} // namespace oomph
 
#endif