Material.h

Go to the documentation of this file.

#ifndef MATERIAL
#define MATERIAL

#include <vector>
#include <string>
#include <iostream>
#include <algorithm>

#include <cmath>
#include <cassert>

#include "Galois/Runtime/PerCPU.h"
class Material {
};

class SimpleMaterial: public Material {
public:

  // (amber) some constants collected here
  static const double EPS = 1.e-6;
  static const double PERT = 1.e-1;
  static const double DET_MIN = 1.e-10;

  static const size_t MAT_SIZE = 9;
  static const size_t NDF = 3;
  static const size_t NDM = 3;

  static const double I_MAT[MAT_SIZE];

  enum ConstRespMode {
    COMPUTE_TANGENTS,
    SKIP_TANGENTS,
  };

  inline SimpleMaterial(double rhoInput = 0) :
    RefRho(rhoInput) {
  }

  virtual ~SimpleMaterial() {
  }

  inline SimpleMaterial(const SimpleMaterial &SM) :
    RefRho(SM.RefRho) {
  }

  virtual SimpleMaterial * clone() const = 0;

  virtual bool getConstitutiveResponse(const std::vector<double>& strain, std::vector<double>& stress, std::vector<double>& tangents
      , const ConstRespMode& mode) const = 0;

  double getDensityInReference() const {
    return RefRho;
  }

  bool getLocalMaterialDensity(const std::vector<double> * strain, double &LocDensity) const {
    assert((*strain).size () == MAT_SIZE);

    // Compute determinant of strain.
    double J = (*strain)[0] * ((*strain)[4] * (*strain)[8] - (*strain)[5] * (*strain)[7]) - (*strain)[1] * ((*strain)[3] * (*strain)[8]
        - (*strain)[5] * (*strain)[6]) + (*strain)[2] * ((*strain)[3] * (*strain)[7] - (*strain)[4] * (*strain)[6]);

    if (fabs(J) < 1.e-10)
      return false;
    else {
      LocDensity = RefRho / J;
      return true;
    }
  }

  virtual const std::string getMaterialName() const = 0;

  static bool consistencyTest(const SimpleMaterial &Smat);

  virtual double getSoundSpeed(void) const = 0;
private:
  double RefRho;
};



class NeoHookean: public SimpleMaterial {

  struct NeoHookenTmpVec {
    static const size_t MAT_SIZE = SimpleMaterial::MAT_SIZE;

    double F[MAT_SIZE];
    double Finv[MAT_SIZE]; 
    double C[MAT_SIZE]; 
    double Cinv[MAT_SIZE]; 
    double S[MAT_SIZE]; 
    double M[MAT_SIZE * MAT_SIZE];

  };

  static GaloisRuntime::PerCPU<NeoHookenTmpVec> perCPUtmpVec;

public:
  NeoHookean(double LambdaInput, double MuInput, double rhoInput = 0) :
    SimpleMaterial(rhoInput), Lambda(LambdaInput), Mu(MuInput) {
  }
  virtual ~NeoHookean() {
  }
  NeoHookean(const NeoHookean &NewMat) :
    SimpleMaterial(NewMat), Lambda(NewMat.Lambda), Mu(NewMat.Mu) {
  }
  virtual NeoHookean * clone() const {
    return new NeoHookean(*this);
  }

  bool getConstitutiveResponse(const std::vector<double>& strain, std::vector<double>& stress, std::vector<double>& tangents
      , const ConstRespMode& mode) const;

  const std::string getMaterialName() const {
    return "NeoHookean";
  }

  double getSoundSpeed(void) const {
    return sqrt((Lambda + 2.0 * Mu) / getDensityInReference());
  }
  ;

private:
  // Lame constants
  double Lambda;
  double Mu;
};

class LinearElasticBase: public SimpleMaterial {
public:
  LinearElasticBase(double irho = 0) :
    SimpleMaterial(irho) {
  }
  virtual ~LinearElasticBase() {
  }
  LinearElasticBase(const LinearElasticBase &NewMat) :
    SimpleMaterial(NewMat) {
  }
  virtual LinearElasticBase * clone() const = 0;

  bool getConstitutiveResponse(const std::vector<double>& strain, std::vector<double>& stress, std::vector<double>& tangents
      , const ConstRespMode& mode) const;

  const std::string getMaterialName() const {
    return "LinearElasticBase";
  }

protected:
  virtual const double getModuli(int i1, int i2, int i3, int i4) const = 0;

};

class IsotropicLinearElastic: public LinearElasticBase {
public:
  static const int DELTA_MAT[][3];

  IsotropicLinearElastic(double iLambda, double imu, double irho = 0) :
    LinearElasticBase(irho), lambda(iLambda), mu(imu) {
  }
  virtual ~IsotropicLinearElastic() {
  }
  IsotropicLinearElastic(const IsotropicLinearElastic &NewMat) :
    LinearElasticBase(NewMat), lambda(NewMat.lambda), mu(NewMat.mu) {
  }
  virtual IsotropicLinearElastic * clone() const {
    return new IsotropicLinearElastic(*this);
  }

  const std::string getMaterialName() const {
    return "IsotropicLinearElastic";
  }

  double getSoundSpeed(void) const {
    return sqrt((lambda + 2.0 * mu) / getDensityInReference());
  }
  ;

protected:
  const double getModuli(int i1, int i2, int i3, int i4) const {
    return lambda * DELTA_MAT[i1][i2] * DELTA_MAT[i3][i4] + mu * (DELTA_MAT[i1][i3] * DELTA_MAT[i2][i4] + DELTA_MAT[i1][i4] * DELTA_MAT[i2][i3]);
  }

private:
  double lambda;
  double mu;
};


#endif

---