young_laplace_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 YoungLaplace elements
#ifndef OOMPH_YOUNGLAPLACE_ELEMENTS_HEADER
#define OOMPH_YOUNGLAPLACE_ELEMENTS_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
 
// OOMPH-LIB headers
#include "../generic/nodes.h"
#include "../generic/Qelements.h"
 
 
namespace oomph
{
  //=============================================================
  //=============================================================
  class YoungLaplaceEquations : public virtual FiniteElement
  {
  public:
    typedef void (*SpineBaseFctPt)(const Vector<double>& x,
                                   Vector<double>& spine_base,
                                   Vector<Vector<double>>& dspine_base);
 
    typedef void (*SpineFctPt)(const Vector<double>& x,
                               Vector<double>& spine,
                               Vector<Vector<double>>& dspine);
 
    YoungLaplaceEquations() : Spine_base_fct_pt(0), Spine_fct_pt(0), Kappa_pt(0)
    {
    }
 
    YoungLaplaceEquations(const YoungLaplaceEquations& dummy) = delete;
 
    void operator=(const YoungLaplaceEquations&) = delete;
 
    virtual inline double u(const unsigned& n) const
    {
      return nodal_value(n, 0);
    }
 
    void output(std::ostream& outfile)
    {
      unsigned n_plot = 5;
      output(outfile, n_plot);
    }
 
    void output(std::ostream& outfile, const unsigned& n_plot);
 
 
    void output_fct(std::ostream& outfile,
                    const unsigned& n_plot,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt);
 
 
    virtual void output_fct(
      std::ostream& outfile,
      const unsigned& n_plot,
      const double& time,
      FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
    {
      throw OomphLibError("These equations are steady => no time dependence",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
 
    void compute_error(std::ostream& outfile,
                       FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
                       double& error,
                       double& norm);
 
    void compute_error(std::ostream& outfile,
                       FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
                       const double& time,
                       double& error,
                       double& norm)
    {
      throw OomphLibError("These equations are steady => no time dependence",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
 
    Data* kappa_pt()
    {
      return Kappa_pt;
    }
 
    bool use_spines() const
    {
      return (Spine_fct_pt != 0);
    }
 
    SpineBaseFctPt& spine_base_fct_pt()
    {
      return Spine_base_fct_pt;
    }
 
 
    SpineFctPt& spine_fct_pt()
    {
      return Spine_fct_pt;
    }
 
    void set_kappa(Data* kappa_pt)
    {
#ifdef PARANOID
      if (kappa_pt->nvalue() != 1)
      {
        throw OomphLibError("kappa must only store a single value!",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Make a local copy
      Kappa_pt = kappa_pt;
 
      // Add to external data and store index in this storage scheme
      Kappa_index = add_external_data(kappa_pt);
    }
 
 
    double get_kappa() const
    {
      if (Kappa_pt == 0)
      {
        return 0.0;
      }
      else
      {
        // Get prescribed kappa value
        return Kappa_pt->value(0);
      }
    }
 
 
    void get_flux(const Vector<double>& s, Vector<double>& flux) const
    {
      // Find out how many nodes there are in the element
      unsigned n_node = nnode();
 
      // Set up memory for the shape (same as test functions)
      Shape psi(n_node);
      DShape dpsidx(n_node, 2);
 
      // Call the derivatives of the shape (same as test functions)
      dshape_eulerian(s, psi, dpsidx);
 
      // Initialise to zero
      for (unsigned j = 0; j < 2; j++)
      {
        flux[j] = 0.0;
      }
 
      // Loop over nodes
      for (unsigned l = 0; l < n_node; l++)
      {
        // Loop over derivative directions
        for (unsigned j = 0; j < 2; j++)
        {
          flux[j] += u(l) * dpsidx(l, j);
        }
      }
    }
 
 
    inline virtual void get_spine_base(
      const Vector<double>& x,
      Vector<double>& spine_base,
      Vector<Vector<double>>& dspine_base) const
    {
      // If no spine function has been set, default to vertical spines
      // emanating from x[0](,x[1])
      if (Spine_base_fct_pt == 0)
      {
        for (unsigned i = 0; i < 3; i++)
        {
          spine_base[i] = x[i];
          for (unsigned j = 0; j < 2; j++)
          {
            dspine_base[i][j] = 0.0;
          }
        }
      }
      else
      {
        // Get spine
        (*Spine_base_fct_pt)(x, spine_base, dspine_base);
      }
    }
 
    inline void get_spine(const Vector<double>& x,
                          Vector<double>& spine,
                          Vector<Vector<double>>& dspine) const
    {
      // If no spine function has been set, default to vertical spines
      // emanating from x[0](,x[1])
      if (Spine_fct_pt == 0)
      {
        // Initialise all to zero
        for (unsigned i = 0; i < 3; i++)
        {
          spine[i] = 0.0;
          for (unsigned j = 0; j < 2; j++)
          {
            dspine[i][j] = 0.0;
          }
        }
        // Overwrite vertical component
        spine[2] = 1.0;
      }
      else
      {
        // Get spine
        (*Spine_fct_pt)(x, spine, dspine);
      }
    }
 
    void position(const Vector<double>& s, Vector<double>& r) const;
 
    void exact_position(const Vector<double>& s,
                        Vector<double>& r,
                        FiniteElement::SteadyExactSolutionFctPt exact_soln_pt);
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals);
 
 
    inline double interpolated_u(const Vector<double>& s) const
    {
      // Find number of nodes
      unsigned n_node = nnode();
 
      // Local shape function
      Shape psi(n_node);
 
      // Find values of shape function
      shape(s, psi);
 
      // Initialise value of u
      double interpolated_u = 0.0;
 
      // Loop over the local nodes and sum
      for (unsigned l = 0; l < n_node; l++)
      {
        interpolated_u += u(l) * psi[l];
      }
 
      return (interpolated_u);
    }
 
 
    unsigned self_test();
 
 
    static void allocate_vector_of_vectors(unsigned n_rows,
                                           unsigned n_cols,
                                           Vector<Vector<double>>& v)
    {
      v.resize(n_rows);
      for (unsigned i = 0; i < n_rows; i++)
      {
        v[i].resize(n_cols);
        for (unsigned j = 0; j < n_cols; j++)
        {
          v[i][j] = 0.0;
        }
      }
    }
 
    static void scalar_times_vector(const double& lambda,
                                    const Vector<double>& v,
                                    Vector<double>& lambda_times_v)
    {
      unsigned n = v.size();
      for (unsigned i = 0; i < n; i++)
      {
        lambda_times_v[i] = lambda * v[i];
      }
    }
 
 
    static double two_norm(const Vector<double>& v)
    {
      double norm = 0.0;
      unsigned n = v.size();
      for (unsigned i = 0; i < n; i++)
      {
        norm += v[i] * v[i];
      }
      return sqrt(norm);
    }
 
 
    static double scalar_product(const Vector<double>& v1,
                                 const Vector<double>& v2)
    {
      double scalar = 0.0;
      unsigned n = v1.size();
      for (unsigned i = 0; i < n; i++)
      {
        scalar += v1[i] * v2[i];
      }
      return scalar;
    }
 
    static void cross_product(const Vector<double>& v1,
                              const Vector<double>& v2,
                              Vector<double>& v_cross)
    {
#ifdef PARANOID
      if ((v1.size() != v2.size()) || (v1.size() != 3))
      {
        throw OomphLibError("Vectors must be of dimension 3 for cross-product!",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      v_cross[0] = v1[1] * v2[2] - v1[2] * v2[1];
      v_cross[1] = v1[2] * v2[0] - v1[0] * v2[2];
      v_cross[2] = v1[0] * v2[1] - v1[1] * v2[0];
    }
 
    static void vector_sum(const Vector<double>& v1,
                           const Vector<double>& v2,
                           Vector<double>& vs)
    {
      unsigned n = v1.size();
      for (unsigned i = 0; i < n; i++)
      {
        vs[i] = v1[i] + v2[i];
      }
    }
 
  protected:
    virtual inline int u_local_eqn(const unsigned& n)
    {
      return nodal_local_eqn(n, 0);
    }
 
    unsigned Kappa_index;
 
    SpineBaseFctPt Spine_base_fct_pt;
 
    SpineFctPt Spine_fct_pt;
 
  private:
    Data* Kappa_pt;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  template<unsigned NNODE_1D>
  class QYoungLaplaceElement : public virtual QElement<2, NNODE_1D>,
                               public virtual YoungLaplaceEquations
  {
  private:
    static const unsigned Initial_Nvalue[];
 
  public:
    QYoungLaplaceElement() : QElement<2, NNODE_1D>(), YoungLaplaceEquations() {}
 
    QYoungLaplaceElement(const QYoungLaplaceElement<NNODE_1D>& dummy) = delete;
 
    void operator=(const QYoungLaplaceElement<NNODE_1D>&) = delete;
 
    inline unsigned required_nvalue(const unsigned& n) const
    {
      return Initial_Nvalue[n];
    }
 
    void output(std::ostream& outfile)
    {
      YoungLaplaceEquations::output(outfile);
    }
 
 
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      YoungLaplaceEquations::output(outfile, n_plot);
    }
 
 
    void output_fct(std::ostream& outfile,
                    const unsigned& n_plot,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
    {
      YoungLaplaceEquations::output_fct(outfile, n_plot, exact_soln_pt);
    }
 
 
    void output_fct(std::ostream& outfile,
                    const unsigned& n_plot,
                    const double& time,
                    FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
    {
      YoungLaplaceEquations::output_fct(outfile, n_plot, time, exact_soln_pt);
    }
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<QYoungLaplaceElement<NNODE_1D>>
    : public virtual QElement<1, NNODE_1D>
  {
  public:
    FaceGeometry() : QElement<1, NNODE_1D>() {}
  };
 
 
 
 
  //========================================================================
  //========================================================================
  class HeightControlElement : public virtual GeneralisedElement
  {
  public:
    HeightControlElement(Node* control_node_pt, double* prescribed_height_pt)
      : GeneralisedElement()
    {
      // Store pointer to prescribed height value
      Prescribed_height_pt = prescribed_height_pt;
 
      // Store pointer to Node at which the height is controlled
      Control_node_pt = control_node_pt;
 
      // Create curvature Data, add it as internal data, and store
      // its number with the internal data storage scheme
      // (i.e. internal_data_pt(Curvature_data_index) will return
      // the pointer to it...).
      Curvature_data_index = add_internal_data(new Data(1));
 
      // Add control_node_pt as external data
      add_external_data(static_cast<Data*>(control_node_pt));
    }
 
    Data*& kappa_pt()
    {
      return internal_data_pt(Curvature_data_index);
    }
 
    void assign_additional_local_eqn_numbers()
    {
      // Get equation number from generic scheme: The height control
      // equation is "the equation for the curvature" which is
      // stored as internal data in the element and has been
      // numbered as such...
      Height_ctrl_local_eqn = internal_local_eqn(Curvature_data_index, 0);
    }
 
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      residuals[Height_ctrl_local_eqn] =
        Control_node_pt->value(0) - (*Prescribed_height_pt);
    }
 
  private:
    double* Prescribed_height_pt;
 
    Node* Control_node_pt;
 
    unsigned Curvature_data_index;
 
    unsigned Height_ctrl_local_eqn;
  };
 
 
} // namespace oomph
 
 
#endif