algebraic_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//====================================================================
#ifndef OOMPH_ALGEBRAIC_ELEMENTS_HEADER
#define OOMPH_ALGEBRAIC_ELEMENTS_HEADER
 
#include "geom_objects.h"
#include "mesh.h"
#include "elements.h"
#include "domain.h"
#include "element_with_moving_nodes.h"
 
namespace oomph
{
  // forward references
  class AlgebraicMesh;
  class AlgebraicElementBase;
  class DummyAlgebraicMesh;
 
 
  // Algebraic nodes
 
 
  //========================================================================
  //========================================================================
  class AlgebraicNode : public Node
  {
  public:
    AlgebraicNode() : Node() {}
 
    AlgebraicNode(const unsigned& n_dim,
                  const unsigned& n_position_type,
                  const unsigned& initial_nvalue)
      : Node(n_dim, n_position_type, initial_nvalue)
    {
#ifdef LEAK_CHECK
      LeakCheckNames::AlgebraicNode_build += 1;
#endif
 
      // Add default node update info
      add_node_update_info(Dummy_node_update_fct_id, // dummy remesh fct ID
                           Dummy_mesh_pt, // dummy mesh
                           Dummy_geom_object_pt, // dummy geom object vector
                           Dummy_ref_value, // dummy  ref vector
                           true); // flag indicating call from
                                  // constructor
    }
 
 
    AlgebraicNode(TimeStepper* 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)
    {
#ifdef LEAK_CHECK
      LeakCheckNames::AlgebraicNode_build += 1;
#endif
 
      // Add default node update info
      add_node_update_info(Dummy_node_update_fct_id, // dummy remesh fct ID
                           Dummy_mesh_pt, // dummy mesh
                           Dummy_geom_object_pt, // dummy geom object vector
                           Dummy_ref_value, // dummy ref vector
                           true); // flag indicating call from
                                  // constructor
    }
 
    virtual ~AlgebraicNode()
    {
#ifdef LEAK_CHECK
      LeakCheckNames::AlgebraicNode_build -= 1;
#endif
    }
 
    AlgebraicNode(const AlgebraicNode&) = delete;
 
    // Commented out broken assignment operator because this can lead to a
    // conflict warning when used in the virtual inheritence hierarchy.
    // Essentially the compiler doesn't realise that two separate
    // implementations of the broken function are the same and so, quite
    // rightly, it shouts.
    /*void operator=(const AlgebraicNode&) = delete;*/
 
 
    void node_update(const bool& update_all_time_levels_for_new_node = false);
 
 
    unsigned nnode_update_fcts()
    {
      // Note: We could read this information out from any one of
      // various maps that that store the information for the
      // different node update functions...
      return Mesh_pt.size();
    }
 
 
    int node_update_fct_id()
    {
      return Default_node_update_fct_id;
    }
 
    void node_update_fct_id(Vector<int>& id)
    {
      // Resize vector
      id.resize(0);
 
      // Loop over all entries and copy them across (again, we could
      // get this information from any of the maps...)
      typedef std::map<int, AlgebraicMesh*>::iterator IT;
      for (IT it = Mesh_pt.begin(); it != Mesh_pt.end(); it++)
      {
        id.push_back(it->first);
      }
    }
 
 
    AlgebraicMesh* mesh_pt()
    {
      return Default_it_mesh_pt->second;
    }
 
 
    AlgebraicMesh* mesh_pt(const int& id)
    {
      return Mesh_pt[id];
    }
 
 
    unsigned ngeom_object(const int& id)
    {
      return Geom_object_pt[id].size();
    }
 
 
    unsigned ngeom_object() const
    {
      return Default_it_geom_object_pt->second.size();
    }
 
 
    Vector<GeomObject*>& vector_geom_object_pt(const int& id)
    {
      return Geom_object_pt[id];
    }
 
 
    Vector<GeomObject*>& vector_geom_object_pt()
    {
      return Default_it_geom_object_pt->second;
    }
 
 
    GeomObject** all_geom_object_pt()
    {
      if (this->ngeom_object() == 0)
      {
        return 0;
      }
      else
      {
        return &(Default_it_geom_object_pt->second[0]);
      }
    }
 
    GeomObject* geom_object_pt(const unsigned& i)
    {
      return Default_it_geom_object_pt->second[i];
    }
 
    unsigned nref_value(const int& id)
    {
      return Ref_value[id].size();
    }
 
 
    unsigned nref_value()
    {
      return Default_it_ref_value->second.size();
    }
 
 
    Vector<double>& vector_ref_value()
    {
      return Default_it_ref_value->second;
    }
 
 
    Vector<double>& vector_ref_value(const int& id)
    {
      return Ref_value[id];
    }
 
 
    double ref_value(const unsigned& i)
    {
      return Default_it_ref_value->second[i];
    }
 
    void add_node_update_info(const int& id,
                              AlgebraicMesh* mesh_pt,
                              const Vector<GeomObject*>& geom_object_pt,
                              const Vector<double>& ref_value,
                              const bool& called_from_constructor = false)
    {
      // Sanity check
      if (id < 0)
      {
        if (!called_from_constructor)
        {
          std::ostringstream error_message;
          error_message << "\nNegative ID, " << id
                        << ", only allowed if called from constructor and\n"
                        << "indicated as such by optional boolean flag."
                        << std::endl;
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // If there's just one entry -- check if it's the default dummy one
      if (Mesh_pt.size() == 1)
      {
        if (Mesh_pt.begin()->second == Dummy_mesh_pt)
        {
          if (Default_it_mesh_pt->second == Dummy_mesh_pt)
          {
            kill_node_update_info(Dummy_node_update_fct_id);
          }
        }
      }
 
      // Now insert the actual info
      Mesh_pt.insert(std::make_pair(id, mesh_pt));
      Geom_object_pt.insert(std::make_pair(id, geom_object_pt));
      Ref_value.insert(std::make_pair(id, ref_value));
 
      // Always use the "first" update fct as default -- can be overwritten
      // outside (usually only done for self test)
      set_default_node_update(Mesh_pt.begin()->first);
    }
 
 
    void add_node_update_info(AlgebraicMesh* mesh_pt,
                              const Vector<GeomObject*>& geom_object_pt,
                              const Vector<double>& ref_value)
    {
      // No update fct id supplied: Use a default assignment of 0.
      unsigned id = 0;
 
      // If there's just one entry -- check if it's the default dummy one
      if (Mesh_pt.size() == 1)
      {
        // Do we still have dummy default assignment stored as the one
        // and only entry?
        if (Mesh_pt.begin()->second == Dummy_mesh_pt)
        {
          if (Default_it_mesh_pt->second == Dummy_mesh_pt)
          {
            kill_node_update_info(Dummy_node_update_fct_id);
          }
        }
      }
 
      // Now insert the actual info
      Mesh_pt.insert(std::make_pair(id, mesh_pt));
      Geom_object_pt.insert(std::make_pair(id, geom_object_pt));
      Ref_value.insert(std::make_pair(id, ref_value));
 
      // Always use the "first" update fct as default -- can be overwritten
      // outside (usually only done for self test)
      set_default_node_update(Mesh_pt.begin()->first);
    }
 
 
    void kill_node_update_info(const int& id = 0)
    {
      Mesh_pt.erase(Mesh_pt.find(id));
      Geom_object_pt.erase(Geom_object_pt.find(id));
      Ref_value.erase(Ref_value.find(id));
    }
 
 
    unsigned self_test();
 
 
  private:
    void set_default_node_update(const int& id)
    {
      // Set default node update fct id
      Default_node_update_fct_id = id;
 
 
      // Set iterators for default entry
 
      // Iterator to default mesh:
      Default_it_mesh_pt = Mesh_pt.find(id);
#ifdef PARANOID
      if (Default_it_mesh_pt == Mesh_pt.end())
      {
        std::ostringstream error_message;
        error_message << "There is no reference mesh for node update fct id"
                      << id << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Iterator to default GeomObject vector
      Default_it_geom_object_pt = Geom_object_pt.find(id);
#ifdef PARANOID
      if (Default_it_geom_object_pt == Geom_object_pt.end())
      {
        std::ostringstream error_message;
        error_message << "There is no Geom_object_pt for node update fct id"
                      << id << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Iterator to default values vector
      Default_it_ref_value = Ref_value.find(id);
#ifdef PARANOID
      if (Default_it_ref_value == Ref_value.end())
      {
        std::ostringstream error_message;
        error_message << "There is no Ref_value for node update fct id" << id
                      << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
 
    std::map<int, AlgebraicMesh*> Mesh_pt;
 
    std::map<int, Vector<GeomObject*>> Geom_object_pt;
 
    std::map<int, Vector<double>> Ref_value;
 
    std::map<int, AlgebraicMesh*>::iterator Default_it_mesh_pt;
 
    std::map<int, Vector<GeomObject*>>::iterator Default_it_geom_object_pt;
 
    std::map<int, Vector<double>>::iterator Default_it_ref_value;
 
    int Default_node_update_fct_id;
 
    static double Max_allowed_difference_between_node_update_fcts;
 
    static int Dummy_node_update_fct_id;
 
    static AlgebraicMesh* Dummy_mesh_pt;
 
    static DummyAlgebraicMesh Dummy_mesh;
 
    static Vector<GeomObject*> Dummy_geom_object_pt;
 
    static Vector<double> Dummy_ref_value;
  };
 
 
  // Algebraic elements
 
 
  //========================================================================
  //========================================================================
  class AlgebraicElementBase
  {
  public:
    AlgebraicElementBase() {}
 
    AlgebraicElementBase(const AlgebraicElementBase&) = delete;
 
    void operator=(const AlgebraicElementBase&) = delete;
 
    void setup_algebraic_node_update(Node*& node_pt,
                                     const Vector<double>& s_father,
                                     FiniteElement* father_el_pt) const;
  };
 
 
  //========================================================================
  //========================================================================
  template<class ELEMENT>
  class AlgebraicElement
    : public ElementWithSpecificMovingNodes<ELEMENT, AlgebraicNode>,
      public AlgebraicElementBase
  {
  public:
    AlgebraicElement()
      : ElementWithSpecificMovingNodes<ELEMENT, AlgebraicNode>(),
        AlgebraicElementBase()
    {
    }
 
    AlgebraicElement(FiniteElement* const& element_pt, const int& face_index)
      : ElementWithSpecificMovingNodes<ELEMENT, AlgebraicNode>(element_pt,
                                                               face_index),
        AlgebraicElementBase()
    {
    }
 
    AlgebraicElement(const AlgebraicElement&) = delete;
 
    /*void operator=(const AlgebraicElement&) = delete;*/
 
 
    ~AlgebraicElement() {}
  };
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class FaceGeometry<AlgebraicElement<ELEMENT>>
    : public virtual FaceGeometry<ELEMENT>
  {
  public:
    FaceGeometry() : FaceGeometry<ELEMENT>() {}
 
  protected:
  };
 
 
  // Algebraic meshes
 
 
  //========================================================================
  //========================================================================
  class AlgebraicMesh : public virtual Mesh
  {
  public:
    AlgebraicMesh()
    {
      add_geom_object_list_pt(0);
    }
 
    AlgebraicMesh(const AlgebraicMesh&) = delete;
 
    /*void operator=(const AlgebraicMesh&) = delete;*/
 
    ~AlgebraicMesh() {}
 
    // Can safely cast the nodes to AlgebraicNodes
    AlgebraicNode* node_pt(const unsigned long& n)
    {
#ifdef PARANOID
      if (!dynamic_cast<AlgebraicNode*>(Node_pt[n]))
      {
        std::ostringstream error_stream;
        error_stream << "Error: Node " << n << "is a "
                     << typeid(Node_pt[n]).name() << ", not an AlgebraicNode"
                     << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Return a cast to the Node_pt
      return (dynamic_cast<AlgebraicNode*>(Node_pt[n]));
    }
 
 
    virtual void algebraic_node_update(const unsigned& t,
                                       AlgebraicNode*& node_pt) = 0;
 
    virtual void update_node_update(AlgebraicNode*& node_pt) = 0;
 
 
    void node_update(const bool& update_all_solid_nodes = false)
    {
#ifdef PARANOID
      if (update_all_solid_nodes)
      {
        std::string error_message =
          "Doesn't make sense to use an AlgebraicMesh with\n";
        error_message +=
          "SolidElements so specifying update_all_solid_nodes=true\n";
        error_message += "doesn't make sense either\n";
 
        throw OomphLibError(
          error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Initial loop over ALL nodes to setup (need to place at least
      // all master nodes before we can update the position of the
      // hanging ones)
      AlgebraicNode* alg_nod_pt = 0;
      unsigned n_node = nnode();
 
      // In parallel there may be no nodes on a particular process
      if (n_node > 0)
      {
        for (unsigned n = 0; n < n_node; n++)
        {
          alg_nod_pt = static_cast<AlgebraicNode*>(node_pt(n));
          alg_nod_pt->node_update();
        }
 
        // Figure out spatial dimension of node
        unsigned n_dim = alg_nod_pt->ndim();
 
        // Now loop over hanging nodes and adjust their nodal positions
        // to reflect the hanging node constraints
        for (unsigned n = 0; n < n_node; n++)
        {
          Node* nod_pt = node_pt(n);
          if (nod_pt->is_hanging())
          {
            // Initialise
            Vector<double> x(n_dim);
            for (unsigned i = 0; i < n_dim; i++)
            {
              x[i] = 0.0;
            }
 
            // Loop over master nodes
            unsigned nmaster = nod_pt->hanging_pt()->nmaster();
            for (unsigned imaster = 0; imaster < nmaster; imaster++)
            {
              // Loop over directions
              for (unsigned i = 0; i < n_dim; i++)
              {
                x[i] += nod_pt->hanging_pt()->master_node_pt(imaster)->x(i) *
                        nod_pt->hanging_pt()->master_weight(imaster);
              }
            }
 
            // Copy across
            for (unsigned i = 0; i < n_dim; i++)
            {
              nod_pt->x(i) = x[i];
            }
            nod_pt->perform_auxiliary_node_update_fct();
          }
        }
      } // end if (n_node>0)
 
#ifdef OOMPH_HAS_MPI
      // Update positions for external halo nodes attached to this mesh
      // Loop over processors
      for (std::map<unsigned, Vector<Node*>>::iterator it =
             External_halo_node_pt.begin();
           it != External_halo_node_pt.end();
           it++)
      {
        int iproc = (*it).first;
        AlgebraicNode* alg_nod_pt = 0;
        unsigned n_ext_halo_node = nexternal_halo_node(iproc);
        // Only act if there are any external halo nodes
        if (n_ext_halo_node > 0)
        {
          for (unsigned n = 0; n < n_ext_halo_node; n++)
          {
            alg_nod_pt =
              static_cast<AlgebraicNode*>(external_halo_node_pt(iproc, n));
            alg_nod_pt->node_update();
          }
 
          // Figure out spatial dimension of node
          unsigned n_dim = alg_nod_pt->ndim();
 
          // Now loop over hanging nodes and adjust their nodal positions
          // to reflect the hanging node constraints
          for (unsigned n = 0; n < n_ext_halo_node; n++)
          {
            Node* nod_pt = external_halo_node_pt(iproc, n);
            if (nod_pt->is_hanging())
            {
              // Initialise
              Vector<double> x(n_dim);
              for (unsigned i = 0; i < n_dim; i++)
              {
                x[i] = 0.0;
              }
 
              // Loop over master nodes
              unsigned nmaster = nod_pt->hanging_pt()->nmaster();
              for (unsigned imaster = 0; imaster < nmaster; imaster++)
              {
                // Loop over directions
                for (unsigned i = 0; i < n_dim; i++)
                {
                  x[i] += nod_pt->hanging_pt()->master_node_pt(imaster)->x(i) *
                          nod_pt->hanging_pt()->master_weight(imaster);
                }
              }
 
              // Copy across
              for (unsigned i = 0; i < n_dim; i++)
              {
                nod_pt->x(i) = x[i];
              }
            }
          }
        }
 
      } // end loop over processors
#endif
    }
 
    unsigned self_test()
    {
      // Initialise
      bool passed = true;
 
      unsigned test = Mesh::self_test();
      if (test != 0)
      {
        passed = false;
      }
 
      // Loop over nodes
      unsigned n_node = nnode();
      for (unsigned n = 0; n < n_node; n++)
      {
        if (static_cast<AlgebraicNode*>(node_pt(n))->self_test() != 0)
        {
          passed = false;
        }
      }
 
      oomph_info << "Done algnode selftest in mesh" << std::endl;
 
      // Return verdict
      if (passed)
      {
        return 0;
      }
      else
      {
        return 1;
      }
    }
 
    void add_geom_object_list_pt(GeomObject* geom_object_pt)
    {
      Geom_object_list_pt.push_back(geom_object_pt);
    }
 
    unsigned ngeom_object_list_pt()
    {
      return Geom_object_list_pt.size();
    }
 
    GeomObject* geom_object_list_pt(const unsigned& i)
    {
      // Probably should be a range check in here...
      return Geom_object_list_pt[i];
    }
 
  private:
    Vector<GeomObject*> Geom_object_list_pt;
  };
 
 
  // Dummy algebraic mesh
 
 
  //========================================================================
  //========================================================================
  class DummyAlgebraicMesh : public virtual AlgebraicMesh
  {
  public:
    DummyAlgebraicMesh() {}
 
    DummyAlgebraicMesh(const DummyAlgebraicMesh&) = delete;
 
    /*void operator=(const DummyAlgebraicMesh&) = delete;*/
 
    virtual void algebraic_node_update(const unsigned& t,
                                       AlgebraicNode*& node_pt)
    {
    }
 
 
    virtual void update_node_update(AlgebraicNode*& node_pt) {}
 
    virtual void setup_algebraic_node_update(AlgebraicNode*& nod_pt) {}
  };
 
 
} // namespace oomph
 
#endif