pml_time_harmonic_elasticity_tensor.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-2025 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 objects representing the fourth-rank elasticity tensor
// for linear elasticity problems
 
// Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_PML_TIME_HARMONIC_ELASTICITY_TENSOR_HEADER
#define OOMPH_PML_TIME_HARMONIC_ELASTICITY_TENSOR_HEADER
 
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include "generic/oomph_utilities.h"
#include "math.h"
#include <complex>
 
namespace oomph
{
  //=====================================================================
  /// A base class that represents the fourth-rank elasticity tensor
  /// \f$E_{ijkl}\f$ defined such that
  /// \f[\tau_{ij} = E_{ijkl} e_{kl},\f]
  /// where \f$e_{ij}\f$ is the infinitessimal (Cauchy) strain tensor
  /// and \f$\tau_{ij}\f$ is the stress tensor. The symmetries of the
  /// tensor are such that
  /// \f[E_{ijkl} = E_{jikl} = E_{ijlk} = E_{klij}\f]
  /// and thus there are relatively few independent components. These
  /// symmetries are included in the definition of the object so that
  /// non-physical symmetries cannot be accidentally imposed.
  //=====================================================================
  class PMLTimeHarmonicElasticityTensor
  {
    /// Translation table from the four indices to the corresponding
    /// independent component
    static const unsigned Index[3][3][3][3];
 
  protected:
    /// Member function that returns the i-th independent component of the
    /// elasticity tensor
    virtual inline std::complex<double> independent_component(
      const unsigned& i) const
    {
      return std::complex<double>(0.0, 0.0);
    }
 
 
    /// Helper range checking function
    /// (Note that this only captures over-runs in 3D but
    /// errors are likely to be caught in evaluation of the
    /// stress and strain tensors anyway...)
    void range_check(const unsigned& i,
                     const unsigned& j,
                     const unsigned& k,
                     const unsigned& l) const
    {
      if ((i > 2) || (j > 2) || (k > 2) || (l > 2))
      {
        std::ostringstream error_message;
        if (i > 2)
        {
          error_message << "Range Error : Index 1 " << i
                        << " is not in the range (0,2)";
        }
        if (j > 2)
        {
          error_message << "Range Error : Index 2 " << j
                        << " is not in the range (0,2)";
        }
 
        if (k > 2)
        {
          error_message << "Range Error : Index 2 " << k
                        << " is not in the range (0,2)";
        }
 
        if (l > 2)
        {
          error_message << "Range Error : Index 4 " << l
                        << " is not in the range (0,2)";
        }
 
        // Throw the error
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
 
 
    /// Empty Constructor
    PMLTimeHarmonicElasticityTensor() {}
 
  public:
    /// Empty virtual Destructor
    virtual ~PMLTimeHarmonicElasticityTensor() {}
 
  public:
    /// Return the appropriate independent component
    /// via the index translation scheme (const version).
    std::complex<double> operator()(const unsigned& i,
                                    const unsigned& j,
                                    const unsigned& k,
                                    const unsigned& l) const
    {
      // Range check
#ifdef PARANOID
      range_check(i, j, k, l);
#endif
      return independent_component(Index[i][j][k][l]);
    }
  };
 
 
  //===================================================================
  /// An isotropic elasticity tensor defined in terms of Young's modulus
  /// and Poisson's ratio. The elasticity tensor is assumed to be
  /// non-dimensionalised on some reference value for Young's modulus
  /// so the value provided to the constructor (if any) is to be
  /// interpreted as the ratio of the actual Young's modulus to the
  /// Young's modulus used to non-dimensionalise the stresses/tractions
  /// in the governing equations.
  //===================================================================
  class PMLTimeHarmonicIsotropicElasticityTensor
    : public PMLTimeHarmonicElasticityTensor
  {
    /// Storage for the independent components of the elasticity tensor
    std::complex<double> C[4];
 
    /// Translation scheme between the 21 independent components of the
    /// general elasticity tensor and the isotropic case
    static const unsigned StaticIndex[21];
 
    /// Poisson's ratio
    double Nu;
 
  public:
    /// Constructor. Passing in the values of the Poisson's ratio
    /// and Young's modulus (interpreted as the ratio of the actual
    /// Young's modulus to the Young's modulus (or other reference stiffness)
    /// used to non-dimensionalise stresses and tractions in the governing
    /// equations).
    PMLTimeHarmonicIsotropicElasticityTensor(const double& nu, const double& E)
      : PMLTimeHarmonicElasticityTensor(), Nu(nu)
    {
      // Set the three indepdent components
      C[0] = std::complex<double>(0.0, 0.0);
      double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
      double mu = E / (2.0 * (1.0 + nu));
      this->set_lame_coefficients(lambda, mu);
    }
 
    /// Constructor. Passing in the value of the Poisson's ratio.
    /// Stresses and tractions in the governing equations are assumed
    /// to have been non-dimensionalised on Young's modulus.
    PMLTimeHarmonicIsotropicElasticityTensor(const double& nu)
      : PMLTimeHarmonicElasticityTensor(), Nu(nu)
    {
      // Set the three indepdent components
      C[0] = std::complex<double>(0.0, 0.0);
 
      // reference value
      double E = 1.0;
      double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
      double mu = E / (2.0 * (1.0 + nu));
      this->set_lame_coefficients(lambda, mu);
    }
 
    /// Poisson's ratio
    double nu() const
    {
      return Nu;
    }
 
    /// Update parameters: Specify values of the Poisson's ratio
    /// and (optionally) Young's modulus (interpreted as the ratio of the actual
    /// Young's modulus to the Young's modulus (or other reference stiffness)
    /// used to non-dimensionalise stresses and tractions in the governing
    /// equations).
    void update_constitutive_parameters(const double& nu, const double& E = 1.0)
    {
      // Set the three indepdent components
      C[0] = std::complex<double>(0.0, 0.0);
      double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
      double mu = E / (2.0 * (1.0 + nu));
      this->set_lame_coefficients(lambda, mu);
    }
 
 
    /// Overload the independent coefficient function
    inline std::complex<double> independent_component(const unsigned& i) const
    {
      return C[StaticIndex[i]];
    }
 
 
  private:
    // Set the values of the lame coefficients
    void set_lame_coefficients(const double& lambda, const double& mu)
    {
      C[1] = lambda + 2.0 * mu;
      C[2] = lambda;
      C[3] = mu;
    }
  };
 
} // namespace oomph
#endif