template_free_poisson.h

//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//           Version 0.90. August 3, 2009.
//LIC//
//LIC// Copyright (C) 2006-2009 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 Poisson elements
#ifndef OOMPH_TF_POISSON_ELEMENTS_HEADER
#define OOMPH_TF_POISSON_ELEMENTS_HEADER


// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif

#include<sstream>

// oomph-lib headers
#include "../../src/generic/Vector.h"
#include "../../src/generic/nodes.h"
#include "../../src/generic/Qelements.h"
#include "../../src/generic/Telements.h"
#include "../../src/generic/oomph_utilities.h"
#include "../../src/generic/oomph_definitions.h"

namespace oomph
{

  //=============================================================
  /// A class for all isoparametric elements that solve the
  /// Poisson equations.
  /// \f[
  /// \frac{\partial^2 u}{\partial x_i^2} = f(x_j)
  /// \f]
  /// This contains the generic maths. Shape functions, geometric
  /// mapping etc. must get implemented in derived class.
  //=============================================================
  class TFPoissonEquations : public virtual FiniteElement
  {

  public:

    /// \short Function pointer to source function fct(x,f(x)) --
    /// x is a Vector!
    typedef void (*PoissonSourceFctPt)(const Vector<double>& x, double& f);


    /// \short Function pointer to gradient of source function  fct(x,g(x)) --
    /// x is a Vector!
    typedef void (*PoissonSourceFctGradientPt)(const Vector<double>& x,
                                               Vector<double>& gradient);


    /// Constructor (must initialise the Source_fct_pt to null)
    TFPoissonEquations() : Source_fct_pt(0), Source_fct_gradient_pt(0)
    {}

    /// Broken copy constructor
    TFPoissonEquations(const TFPoissonEquations& dummy)
    {
      BrokenCopy::broken_copy("TFPoissonEquations");
    }

    /// Broken assignment operator
    void operator=(const TFPoissonEquations&)
    {
      BrokenCopy::broken_assign("TFPoissonEquations");
    }

    /// \short Requires storage for one value.
    unsigned required_nvalue(const unsigned &n) const {return 1;}

    /// \short Return the index at which the unknown value
    /// is stored. The default value, 0, is appropriate for single-physics
    /// problems, when there is only one variable, the value that satisfies
    /// the poisson equation.
    /// In derived multi-physics elements, this function should be overloaded
    /// to reflect the chosen storage scheme. Note that these equations require
    /// that the unknown is always stored at the same index at each node.
    virtual inline unsigned u_index_poisson() const {return 0;}

    /// Output with default number of plot points
    void output(std::ostream &outfile)
    {
      const unsigned n_plot=5;
      output(outfile,n_plot);
    }

    /// \short Output FE representation of soln: x,y,u or x,y,z,u at
    /// n_plot^DIM plot points
    void output(std::ostream &outfile, const unsigned &n_plot);

    /// C_style output with default number of plot points
    void output(FILE* file_pt)
    {
      const unsigned n_plot=5;
      output(file_pt,n_plot);
    }

    /// \short C-style output FE representation of soln: x,y,u or x,y,z,u at
    /// n_plot^DIM plot points
    void output(FILE* file_pt, const unsigned &n_plot);

    /// Output exact soln: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
    void output_fct(std::ostream &outfile, const unsigned &n_plot,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt);

    /// \short Output exact soln: x,y,u_exact or x,y,z,u_exact at
    /// n_plot^DIM plot points (dummy time-dependent version to
    /// keep intel compiler happy)
    virtual void output_fct(std::ostream &outfile, const unsigned &n_plot,
                            const double& time,
                            FiniteElement::UnsteadyExactSolutionFctPt
                            exact_soln_pt)
    {
      throw OomphLibError(
                          "There is no time-dependent output_fct() for Poisson elements ",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }


    /// Get error against and norm of exact solution
    void compute_error(std::ostream &outfile,
                       FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
                       double& error, double& norm);


    /// Dummy, time dependent error checker
    void compute_error(std::ostream &outfile,
                       FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
                       const double& time, double& error, double& norm)
    {
      throw OomphLibError("There is no time-dependent compute_error() for Poisson elements",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }

    /// Access function: Pointer to source function
    PoissonSourceFctPt& source_fct_pt() {return Source_fct_pt;}

    /// Access function: Pointer to source function. Const version
    PoissonSourceFctPt source_fct_pt() const {return Source_fct_pt;}

    /// Access function: Pointer to gradient of source function
    PoissonSourceFctGradientPt& source_fct_gradient_pt()
    {return Source_fct_gradient_pt;}

    /// Access function: Pointer to gradient source function. Const version
    PoissonSourceFctGradientPt source_fct_gradient_pt() const
    {return Source_fct_gradient_pt;}


    /// Get source term at (Eulerian) position x. This function is
    /// virtual to allow overloading in multi-physics problems where
    /// the strength of the source function might be determined by
    /// another system of equations.
    inline virtual void get_source_poisson(const unsigned& ipt,
                                           const Vector<double>& x,
                                           double& source) const
    {
      //If no source function has been set, return zero
      if(Source_fct_pt==0) {source = 0.0;}
      else
        {
          // Get source strength
          (*Source_fct_pt)(x,source);
        }
    }


    /// Get gradient of source term at (Eulerian) position x. This function is
    /// virtual to allow overloading in multi-physics problems where
    /// the strength of the source function might be determined by
    /// another system of equations. Computed via function pointer
    /// (if set) or by finite differencing (default)
    inline virtual void get_source_gradient_poisson(
                                                    const unsigned& ipt,
                                                    const Vector<double>& x,
                                                    Vector<double>& gradient) const
    {
      //If no gradient function has been set, FD it
      if(Source_fct_gradient_pt==0)
        {
          // Reference value
          double source=0.0;
          get_source_poisson(ipt,x,source);

          // FD it
          double eps_fd=GeneralisedElement::Default_fd_jacobian_step;
          double source_pls=0.0;
          Vector<double> x_pls(x);
          for (unsigned i=0;i<nodal_dimension();i++)
            {
              x_pls[i]+=eps_fd;
              get_source_poisson(ipt,x_pls,source_pls);
              gradient[i]=(source_pls-source)/eps_fd;
              x_pls[i]=x[i];
            }
        }
      else
        {
          // Get gradient
          (*Source_fct_gradient_pt)(x,gradient);
        }
    }




    /// Get flux: flux[i] = du/dx_i
    void get_flux(const Vector<double>& s, Vector<double>& flux) const
    {
      //Find out how many nodes there are in the element
      const unsigned n_node = nnode();

      //Get the index at which the unknown is stored
      const unsigned u_nodal_index = u_index_poisson();

      //Set up memory for the shape and test functions
      Shape psi(n_node);
      DShape dpsidx(n_node,nodal_dimension());

      //Call the derivatives of the shape and test functions
      dshape_eulerian(s,psi,dpsidx);

      //Initialise to zero
      for(unsigned j=0;j<nodal_dimension();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<nodal_dimension();j++)
            {
              flux[j] += this->nodal_value(l,u_nodal_index)*dpsidx(l,j);
            }
        }
    }


    /// Add the element's contribution to its residual vector (wrapper)
    void fill_in_contribution_to_residuals(Vector<double> &residuals)
    {
      //Call the generic residuals function with flag set to 0
      //using a dummy matrix argument
      fill_in_generic_residual_contribution_poisson(
                                                    residuals,GeneralisedElement::Dummy_matrix,0);
    }


    /// Add the element's contribution to its residual vector and
    /// element Jacobian matrix (wrapper)
    void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                          DenseMatrix<double> &jacobian)
    {
      //Call the generic routine with the flag set to 1
      fill_in_generic_residual_contribution_poisson(residuals,jacobian,1);
    }



    /// \short Return FE representation of function value u_poisson(s)
    /// at local coordinate s
    inline double interpolated_u_poisson(const Vector<double> &s) const
    {
      //Find number of nodes
      const unsigned n_node = nnode();

      //Get the index at which the poisson unknown is stored
      const unsigned u_nodal_index = u_index_poisson();

      //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 += this->nodal_value(l,u_nodal_index)*psi[l];
        }

      return(interpolated_u);
    }


    /// \short Compute derivatives of elemental residual vector with respect
    /// to nodal coordinates. Overwrites default implementation in
    /// FiniteElement base class.
    /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
    virtual void get_dresidual_dnodal_coordinates(RankThreeTensor<double>&
                                                  dresidual_dnodal_coordinates);

    /// \short Self-test: Return 0 for OK
    unsigned self_test();


  protected:

    /// \short Shape/test functions and derivs w.r.t. to global coords at
    /// local coord. s; return  Jacobian of mapping
    virtual double dshape_and_dtest_eulerian_poisson(const Vector<double> &s,
                                                     Shape &psi,
                                                     DShape &dpsidx, Shape &test,
                                                     DShape &dtestdx) const=0;


    /// \short Shape/test functions and derivs w.r.t. to global coords at
    /// integration point ipt; return  Jacobian of mapping
    virtual double dshape_and_dtest_eulerian_at_knot_poisson(const unsigned &ipt,
                                                             Shape &psi,
                                                             DShape &dpsidx,
                                                             Shape &test,
                                                             DShape &dtestdx)
      const=0;

    /// \short Shape/test functions and derivs w.r.t. to global coords at
    /// integration point ipt; return Jacobian of mapping (J). Also compute
    /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
    virtual double dshape_and_dtest_eulerian_at_knot_poisson(
                                                             const unsigned &ipt,
                                                             Shape &psi,
                                                             DShape &dpsidx,
                                                             RankFourTensor<double> &d_dpsidx_dX,
                                                             Shape &test,
                                                             DShape &dtestdx,
                                                             RankFourTensor<double> &d_dtestdx_dX,
                                                             DenseMatrix<double> &djacobian_dX) const=0;

    /// \short Compute element residual Vector only (if flag=and/or element
    /// Jacobian matrix
    virtual void fill_in_generic_residual_contribution_poisson(
                                                               Vector<double> &residuals, DenseMatrix<double> &jacobian,
                                                               const unsigned& flag);

    /// Pointer to source function:
    PoissonSourceFctPt Source_fct_pt;

    /// Pointer to gradient of source function
    PoissonSourceFctGradientPt Source_fct_gradient_pt;

  };






  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////



  //======================================================================
  /// QTFPoissonElement elements are linear/quadrilateral/brick-shaped
  /// Poisson elements with isoparametric interpolation for the function.
  //======================================================================
  template <unsigned DIM, unsigned NNODE_1D>
  class QTFPoissonElement : public virtual QElement<DIM,NNODE_1D>,
                          public virtual TFPoissonEquations
  {

  public:


    ///\short  Constructor: Call constructors for QElement and
    /// Poisson equations
    QTFPoissonElement() : QElement<DIM,NNODE_1D>(), TFPoissonEquations() {}

    /// Broken copy constructor
    QTFPoissonElement(const QTFPoissonElement<DIM,NNODE_1D>& dummy)
    {
      BrokenCopy::broken_copy("QTFPoissonElement");
    }

    /// Broken assignment operator
    void operator=(const QTFPoissonElement<DIM,NNODE_1D>&)
    {
      BrokenCopy::broken_assign("QTFPoissonElement");
    }


    /// \short  Required  # of `values' (pinned or dofs)
    /// at node n
    inline unsigned required_nvalue(const unsigned &n) const
    {return 1;}

    /// \short Output function:
    ///  x,y,u   or    x,y,z,u
    void output(std::ostream &outfile)
    {TFPoissonEquations::output(outfile);}


    ///  \short Output function:
    ///   x,y,u   or    x,y,z,u at n_plot^DIM plot points
    void output(std::ostream &outfile, const unsigned &n_plot)
    {TFPoissonEquations::output(outfile,n_plot);}


    /// \short C-style output function:
    ///  x,y,u   or    x,y,z,u
    void output(FILE* file_pt)
    {TFPoissonEquations::output(file_pt);}


    ///  \short C-style output function:
    ///   x,y,u   or    x,y,z,u at n_plot^DIM plot points
    void output(FILE* file_pt, const unsigned &n_plot)
    {TFPoissonEquations::output(file_pt,n_plot);}


    /// \short Output function for an exact solution:
    ///  x,y,u_exact   or    x,y,z,u_exact at n_plot^DIM plot points
    void output_fct(std::ostream &outfile, const unsigned &n_plot,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
    {TFPoissonEquations::output_fct(outfile,n_plot,exact_soln_pt);}



    /// \short Output function for a time-dependent exact solution.
    ///  x,y,u_exact   or    x,y,z,u_exact at n_plot^DIM plot points
    /// (Calls the steady version)
    void output_fct(std::ostream &outfile, const unsigned &n_plot,
                    const double& time,
                    FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
    {TFPoissonEquations::output_fct(outfile,n_plot,time,exact_soln_pt);}


  protected:

    /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
    inline double dshape_and_dtest_eulerian_poisson(
                                                    const Vector<double> &s, Shape &psi, DShape &dpsidx,
                                                    Shape &test, DShape &dtestdx) const;


    /// \short Shape, test functions & derivs. w.r.t. to global coords. at
    /// integration point ipt. Return Jacobian.
    inline double dshape_and_dtest_eulerian_at_knot_poisson(const unsigned& ipt,
                                                            Shape &psi,
                                                            DShape &dpsidx,
                                                            Shape &test,
                                                            DShape &dtestdx)
      const;

    /// \short Shape/test functions and derivs w.r.t. to global coords at
    /// integration point ipt; return Jacobian of mapping (J). Also compute
    /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
    inline double dshape_and_dtest_eulerian_at_knot_poisson(
                                                            const unsigned &ipt,
                                                            Shape &psi,
                                                            DShape &dpsidx,
                                                            RankFourTensor<double> &d_dpsidx_dX,
                                                            Shape &test,
                                                            DShape &dtestdx,
                                                            RankFourTensor<double> &d_dtestdx_dX,
                                                            DenseMatrix<double> &djacobian_dX) const;

  };




  //Inline functions:


  //======================================================================
  /// Define the shape functions and test functions and derivatives
  /// w.r.t. global coordinates and return Jacobian of mapping.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double QTFPoissonElement<DIM,NNODE_1D>::
  dshape_and_dtest_eulerian_poisson(const Vector<double> &s,
                                    Shape &psi,
                                    DShape &dpsidx,
                                    Shape &test,
                                    DShape &dtestdx) const
  {
    //Call the geometrical shape functions and derivatives
    const double J = this->dshape_eulerian(s,psi,dpsidx);

    //Set the test functions equal to the shape functions
    test = psi;
    dtestdx= dpsidx;

    //Return the jacobian
    return J;
  }




  //======================================================================
  /// Define the shape functions and test functions and derivatives
  /// w.r.t. global coordinates and return Jacobian of mapping.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double QTFPoissonElement<DIM,NNODE_1D>::
  dshape_and_dtest_eulerian_at_knot_poisson(
                                            const unsigned &ipt,
                                            Shape &psi,
                                            DShape &dpsidx,
                                            Shape &test,
                                            DShape &dtestdx) const
  {
    //Call the geometrical shape functions and derivatives
    const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);

    //Set the pointers of the test functions
    test = psi;
    dtestdx = dpsidx;

    //Return the jacobian
    return J;
  }



  //======================================================================
  /// Define the shape functions (psi) and test functions (test) and
  /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
  /// and return Jacobian of mapping (J). Additionally compute the
  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double QTFPoissonElement<DIM,NNODE_1D>::
  dshape_and_dtest_eulerian_at_knot_poisson(
                                            const unsigned &ipt,
                                            Shape &psi,
                                            DShape &dpsidx,
                                            RankFourTensor<double> &d_dpsidx_dX,
                                            Shape &test,
                                            DShape &dtestdx,
                                            RankFourTensor<double> &d_dtestdx_dX,
                                            DenseMatrix<double> &djacobian_dX) const
  {
    // Call the geometrical shape functions and derivatives
    const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
                                                   djacobian_dX,d_dpsidx_dX);

    // Set the pointers of the test functions
    test = psi;
    dtestdx = dpsidx;
    d_dtestdx_dX = d_dpsidx_dX;

    //Return the jacobian
    return J;
  }



  //======================================================================
  /// TTFPoissonElement elements are linear/quadrilateral/brick-shaped
  /// Poisson elements with isoparametric interpolation for the function.
  //======================================================================
  template <unsigned DIM, unsigned NNODE_1D>
  class TTFPoissonElement : public virtual TElement<DIM,NNODE_1D>,
                            public virtual TFPoissonEquations
  {

  public:


    ///\short  Constructor: Call constructors for TElement and
    /// Poisson equations
    TTFPoissonElement() : TElement<DIM,NNODE_1D>(), TFPoissonEquations()
    {}

    /// Broken copy constructor
    TTFPoissonElement(const TTFPoissonElement<DIM,NNODE_1D>& dummy)
    {
      BrokenCopy::broken_copy("TTFPoissonElement");
    }

    /// Broken assignment operator
    void operator=(const TTFPoissonElement<DIM,NNODE_1D>&)
    {
      BrokenCopy::broken_assign("TTFPoissonElement");
    }


    /// \short  Required  # of `values' (pinned or dofs)
    /// at node n
    inline unsigned required_nvalue(const unsigned &n) const
    {return 1;}

    /// \short Output function:
    ///  x,y,u   or    x,y,z,u
    void output(std::ostream &outfile)
    {TFPoissonEquations::output(outfile);}


    ///  \short Output function:
    ///   x,y,u   or    x,y,z,u at n_plot^DIM plot points
    void output(std::ostream &outfile, const unsigned &n_plot)
    {TFPoissonEquations::output(outfile,n_plot);}


    /// \short C-style output function:
    ///  x,y,u   or    x,y,z,u
    void output(FILE* file_pt)
    {TFPoissonEquations::output(file_pt);}


    ///  \short C-style output function:
    ///   x,y,u   or    x,y,z,u at n_plot^DIM plot points
    void output(FILE* file_pt, const unsigned &n_plot)
    {TFPoissonEquations::output(file_pt,n_plot);}


    /// \short Output function for an exact solution:
    ///  x,y,u_exact   or    x,y,z,u_exact at n_plot^DIM plot points
    void output_fct(std::ostream &outfile, const unsigned &n_plot,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
    {TFPoissonEquations::output_fct(outfile,n_plot,exact_soln_pt);}



    /// \short Output function for a time-dependent exact solution.
    ///  x,y,u_exact   or    x,y,z,u_exact at n_plot^DIM plot points
    /// (Calls the steady version)
    void output_fct(std::ostream &outfile, const unsigned &n_plot,
                    const double& time,
                    FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
    {TFPoissonEquations::output_fct(outfile,n_plot,time,exact_soln_pt);}


  protected:

    /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
    inline double dshape_and_dtest_eulerian_poisson(
                                                    const Vector<double> &s, Shape &psi, DShape &dpsidx,
                                                    Shape &test, DShape &dtestdx) const;


    /// \short Shape, test functions & derivs. w.r.t. to global coords. at
    /// integration point ipt. Return Jacobian.
    inline double dshape_and_dtest_eulerian_at_knot_poisson(const unsigned& ipt,
                                                            Shape &psi,
                                                            DShape &dpsidx,
                                                            Shape &test,
                                                            DShape &dtestdx)
      const;

    /// \short Shape/test functions and derivs w.r.t. to global coords at
    /// integration point ipt; return Jacobian of mapping (J). Also compute
    /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
    inline double dshape_and_dtest_eulerian_at_knot_poisson(
                                                            const unsigned &ipt,
                                                            Shape &psi,
                                                            DShape &dpsidx,
                                                            RankFourTensor<double> &d_dpsidx_dX,
                                                            Shape &test,
                                                            DShape &dtestdx,
                                                            RankFourTensor<double> &d_dtestdx_dX,
                                                            DenseMatrix<double> &djacobian_dX) const;

  };




  //Inline functions:


  //======================================================================
  /// Define the shape functions and test functions and derivatives
  /// w.r.t. global coordinates and return Jacobian of mapping.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double TTFPoissonElement<DIM,NNODE_1D>::
  dshape_and_dtest_eulerian_poisson(const Vector<double> &s,
                                    Shape &psi,
                                    DShape &dpsidx,
                                    Shape &test,
                                    DShape &dtestdx) const
  {
    //Call the geometrical shape functions and derivatives
    const double J = this->dshape_eulerian(s,psi,dpsidx);

    //Set the test functions equal to the shape functions
    test = psi;
    dtestdx= dpsidx;

    //Return the jacobian
    return J;
  }




  //======================================================================
  /// Define the shape functions and test functions and derivatives
  /// w.r.t. global coordinates and return Jacobian of mapping.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double TTFPoissonElement<DIM,NNODE_1D>::
  dshape_and_dtest_eulerian_at_knot_poisson(
                                            const unsigned &ipt,
                                            Shape &psi,
                                            DShape &dpsidx,
                                            Shape &test,
                                            DShape &dtestdx) const
  {
    //Call the geometrical shape functions and derivatives
    const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);

    //Set the pointers of the test functions
    test = psi;
    dtestdx = dpsidx;

    //Return the jacobian
    return J;
  }



  //======================================================================
  /// Define the shape functions (psi) and test functions (test) and
  /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
  /// and return Jacobian of mapping (J). Additionally compute the
  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double TTFPoissonElement<DIM,NNODE_1D>::
  dshape_and_dtest_eulerian_at_knot_poisson(
                                            const unsigned &ipt,
                                            Shape &psi,
                                            DShape &dpsidx,
                                            RankFourTensor<double> &d_dpsidx_dX,
                                            Shape &test,
                                            DShape &dtestdx,
                                            RankFourTensor<double> &d_dtestdx_dX,
                                            DenseMatrix<double> &djacobian_dX) const
  {
    // Call the geometrical shape functions and derivatives
    const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
                                                   djacobian_dX,d_dpsidx_dX);

    // Set the pointers of the test functions
    test = psi;
    dtestdx = dpsidx;
    d_dtestdx_dX = d_dpsidx_dX;

    //Return the jacobian
    return J;
  }





}

#endif