text: // 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-2021 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_TIME_HARMONIC_ELASTICITY_TENSOR_HEADER
  #define OOMPH_TIME_HARMONIC_ELASTICITY_TENSOR_HEADER
  
  // Config header generated by autoconfig
  #ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
  #endif
  
  #include "../generic/oomph_utilities.h"
  
  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 TimeHarmonicElasticityTensor
    {
      /// 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 double independent_component(const unsigned& i) const
      {
        return 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
      TimeHarmonicElasticityTensor() {}
  
    public:
      /// Empty virtual Destructor
      virtual ~TimeHarmonicElasticityTensor() {}
  
    public:
      /// Return the appropriate independent component
      /// via the index translation scheme (const version).
      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 TimeHarmonicIsotropicElasticityTensor
      : public TimeHarmonicElasticityTensor
    {
      // Storage for the independent components of the elasticity tensor
      double C[4];
  
      // Translation scheme between the 21 independent components of the general
      // elasticity tensor and the isotropic case
      static const unsigned StaticIndex[21];
  
    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).
      TimeHarmonicIsotropicElasticityTensor(const double& nu, const double& E)
        : TimeHarmonicElasticityTensor()
      {
        // Set the three indepdent components
        C[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.
      TimeHarmonicIsotropicElasticityTensor(const double& nu)
        : TimeHarmonicElasticityTensor()
      {
        // Set the three indepdent components
        C[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);
      }
  
      /// 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] = 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 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
  