materialDeck.h

Go to the documentation of this file.

/*
 * -------------------------------------------
 * Copyright (c) 2021 - 2026 Prashant K. Jha
 * -------------------------------------------
 * PeriDEM https://github.com/prashjha/PeriDEM
 *
 * Distributed under the Boost Software License, Version 1.0. (See accompanying
 * file LICENSE)
 */
 
#ifndef INP_MATERIALDECK_H
#define INP_MATERIALDECK_H
 
#include "deckField.h"
#include "util/io.h"
#include "util/json.h"
#include <cmath>
#include <string>
#include <vector>
 
namespace inp {
 
struct MatData {
 
  double d_E;
 
  double d_G;
 
  double d_K;
 
  double d_nu;
 
  double d_lambda;
 
  double d_mu;
 
  double d_KIc;
 
  double d_Gc;
 
  MatData()
      : d_E(-1.), d_G(-1.), d_K(-1.), d_nu(-1.), d_lambda(-1.), d_mu(-1.),
        d_KIc(-1.), d_Gc(-1.){};
 
  MatData(const MatData &md)
      : d_E(md.d_E), d_G(md.d_G), d_K(md.d_K), d_nu(md.d_nu),
        d_lambda(md.d_lambda), d_mu(md.d_mu), d_KIc(md.d_KIc), d_Gc(md.d_Gc){};
 
  static const std::vector<Field<MatData>> &fields() {
    static const std::vector<Field<MatData>> f = {
        field(&MatData::d_E, "E", -1., "Young's modulus"),
        field(&MatData::d_G, "G", -1., "Shear modulus"),
        field(&MatData::d_K, "K", -1., "Bulk modulus"),
        field(&MatData::d_lambda, "Lambda", -1., "First Lame parameter"),
        field(&MatData::d_mu, "Mu", -1., "Second Lame parameter"),
        field(&MatData::d_nu, "Poisson_Ratio", -1., "Poisson ratio"),
        field(&MatData::d_Gc, "Gc", -1., "Critical energy release rate"),
        field(&MatData::d_KIc, "KIc", -1., "Critical stress intensity factor"),
    };
    return f;
  }
 
  void readFromJson(const json &j) {
    if (j.empty())
      return;
    readFields(*this, j, fields());
  }
 
  std::string printStr(int nt = 0, int lvl = 0) const {
 
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "------- MatData --------" << std::endl << std::endl;
    oss << tabS << "Young's modulus = " << d_E << std::endl;
    oss << tabS << "Shear modulus = " << d_G << std::endl;
    oss << tabS << "Bulk modulus = " << d_K << std::endl;
    oss << tabS << "Poisson ratio = " << d_nu << std::endl;
    oss << tabS << "Lame parameter Lambda = " << d_lambda << std::endl;
    oss << tabS << "Lame parameter Mu = " << d_mu << std::endl;
    oss << tabS << "Critical stress intensity factor = " << d_KIc << std::endl;
    oss << tabS << "Critical energy release rate = " << d_Gc << std::endl;
    oss << tabS << std::endl;
 
    return oss.str();
  }
 
  void print(int nt = 0, int lvl = 0) const { std::cout << printStr(nt, lvl); }
 
  double toNu(double lambda, double mu) { return lambda * 0.5 / (lambda + mu); }
 
  double toNuEG(double E, double G) { return E * 0.5 / G - 1.; }
 
  double toE(double K, double nu) { return K * (3. * (1. - 2. * nu)); }
 
  double toK(double E, double nu) { return E / (3. * (1. - 2. * nu)); }
 
  double toLambdaE(double E, double nu) {
    return E * nu / ((1. + nu) * (1. - 2. * nu));
  }
 
  double toLambdaK(double K, double nu) { return 3. * K * nu / (1. + nu); }
 
  double toGE(double E, double nu) { return E / (2. * (1. + nu)); }
 
  double toGK(double K, double nu) {
    return 3. * K * (1. - 2. * nu) / (2. * (1. + nu));
  }
 
  double toELambda(double lambda, double nu) {
    return lambda * (1. + nu) * (1. - 2. * nu) / nu;
  }
 
  double toGc(double KIc, double nu, double E) { return KIc * KIc / E; }
 
  double toKIc(double Gc, double nu, double E) { return std::sqrt(Gc * E); }
 
};
 
struct MaterialDeck {
 
  bool d_isPlaneStrain;
 
  std::string d_materialType;
 
  size_t d_bondPotentialType;
 
  size_t d_statePotentialType;
 
  size_t d_influenceFnType;
 
  std::vector<double> d_bondPotentialParams;
 
  std::vector<double> d_statePotentialParams;
 
  std::vector<double> d_influenceFnParams;
 
  bool d_irreversibleBondBreak;
 
  bool d_stateContributionFromBrokenBond;
 
  double d_checkScFactor;
 
  bool d_computeParamsFromElastic;
 
  inp::MatData d_matData;
 
  double d_density;
 
  double d_horizon;
 
  double d_horizonMeshRatio;
 
  MaterialDeck(const json &j = json({}))
      : d_materialType(""),
        d_isPlaneStrain(false), d_bondPotentialType(0), d_statePotentialType(0),
        d_influenceFnType(0), d_irreversibleBondBreak(true),
        d_stateContributionFromBrokenBond(true), d_checkScFactor(1.),
        d_computeParamsFromElastic(true), d_matData(inp::MatData()),
        d_density(1.), d_horizon(-1.), d_horizonMeshRatio(-1.) {
    readFromJson(j);
  };
 
  MaterialDeck(std::string materialType, bool isPlainStrain = false,
               double horizon = -1., double horizonMeshRatio = -1., double density = 1.,
               double K = 0., double G = 0., double Gc = 0., bool computeParamsFromElastic = true,
               size_t influenceFnType = 0)
      : d_materialType(materialType),
        d_isPlaneStrain(isPlainStrain), d_bondPotentialType(0), d_statePotentialType(0),
        d_influenceFnType(influenceFnType), d_irreversibleBondBreak(true),
        d_stateContributionFromBrokenBond(true), d_checkScFactor(1.),
        d_computeParamsFromElastic(computeParamsFromElastic), d_matData(inp::MatData()),
        d_density(density), d_horizon(horizon), d_horizonMeshRatio(horizonMeshRatio) {
 
    d_matData.d_K = K;
    d_matData.d_G = G;
    d_matData.d_Gc = Gc;
  };
 
  MaterialDeck(const MaterialDeck &md)
      : d_materialType(md.d_materialType),
        d_isPlaneStrain(md.d_isPlaneStrain),
        d_bondPotentialType(md.d_bondPotentialType),
        d_statePotentialType(md.d_statePotentialType),
        d_influenceFnType(md.d_influenceFnType),
        d_bondPotentialParams(md.d_bondPotentialParams),
        d_statePotentialParams(md.d_statePotentialParams),
        d_influenceFnParams(md.d_influenceFnParams),
        d_irreversibleBondBreak(md.d_irreversibleBondBreak),
        d_stateContributionFromBrokenBond(md.d_stateContributionFromBrokenBond),
        d_checkScFactor(md.d_checkScFactor),
        d_computeParamsFromElastic(md.d_computeParamsFromElastic),
        d_matData(md.d_matData), d_density(md.d_density),
        d_horizon(md.d_horizon), d_horizonMeshRatio(md.d_horizonMeshRatio){};
 
  static const std::vector<Field<MaterialDeck>> &fields() {
    static const std::vector<Field<MaterialDeck>> f = {
        field(&MaterialDeck::d_materialType, "Type", std::string(),
              "Peridynamic material model"),
        field(&MaterialDeck::d_density, "Density", 1., "Mass density"),
        field(&MaterialDeck::d_computeParamsFromElastic,
              "Compute_From_Classical", true,
              "Derive the bond parameters from the elastic constants"),
        field<MaterialDeck, bool>(&MaterialDeck::d_isPlaneStrain,
                                  "Is_Plane_Strain", false,
                                  "Plane strain rather than plane stress", {},
                                  [](const bool &v) { return v; }),
        field<MaterialDeck, double>(&MaterialDeck::d_horizon, "Horizon", -1.,
                                    "Peridynamic horizon", {},
                                    [](const double &v) { return v > 0.; }),
        field<MaterialDeck, double>(&MaterialDeck::d_horizonMeshRatio,
                                    "Horizon_Mesh_Ratio", -1.,
                                    "Horizon as a multiple of the mesh size",
                                    {},
                                    [](const double &v) { return v > 0.; }),
        field<MaterialDeck, std::vector<double>>(
            &MaterialDeck::d_bondPotentialParams, "Bond_Potential_Params", {},
            "Bond parameters given directly rather than derived", {},
            [](const std::vector<double> &v) { return !v.empty(); }),
        field<MaterialDeck, std::vector<double>>(
            &MaterialDeck::d_statePotentialParams, "State_Potential_Params",
            {}, "State parameters given directly rather than derived", {},
            [](const std::vector<double> &v) { return !v.empty(); }),
    };
    return f;
  }
 
  static const std::vector<Field<MaterialDeck>> &influenceFields() {
    static const std::vector<Field<MaterialDeck>> f = {
        field(&MaterialDeck::d_influenceFnType, "Type", size_t(0),
              "Influence function used in the bond force"),
        field<MaterialDeck, std::vector<double>>(
            &MaterialDeck::d_influenceFnParams, "Parameters", {},
            "Parameters of the influence function", {},
            [](const std::vector<double> &v) { return !v.empty(); }),
    };
    return f;
  }
 
  static const std::vector<OneOf> &groups() {
    static const std::vector<OneOf> g = {
        {"the horizon", {"Horizon", "Horizon_Mesh_Ratio"}},
    };
    return g;
  }
 
  static json getExampleJson(const json &given = json::object()) {
    const std::vector<std::string> extra = {"E", "K", "G", "Gc",
                                            "Influence_Function"};
    json j = applyGiven(given, fields(), extra);
    checkGroups(j, groups());
    return j;
  }
 
  static json getExampleJson(const std::string &) = delete;
  static json getExampleJson(const char *) = delete;
 
  void readFromJson(const json &j) {
    if (j.empty())
      return;
 
    // read mat data
    d_matData.readFromJson(j);
 
    readFields(*this, j, fields());
 
    // Older decks spell this Is_Plain_Strain. The table reads the current
    // spelling, so the older one is applied when it alone is present.
    if (j.find("Is_Plane_Strain") == j.end() &&
        j.find("Is_Plain_Strain") != j.end())
      d_isPlaneStrain = j.at("Is_Plain_Strain").get<bool>();
 
    checkGroups(j, groups());
 
    if (j.find("Influence_Function") != j.end())
      readFields(*this, j.at("Influence_Function"), influenceFields());
    if (!d_computeParamsFromElastic && d_bondPotentialParams.empty())
      throw std::runtime_error(
          "Compute_From_Classical is false but Bond_Potential_Params is missing.");
  }
 
  std::string printStr(int nt = 0, int lvl = 0) const {
 
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "------- MaterialDeck --------" << std::endl << std::endl;
    oss << tabS << "Is plain strain = " << d_isPlaneStrain << std::endl;
    oss << tabS << "Material type = " << d_materialType << std::endl;
    oss << tabS << "Bond potential type = " << d_bondPotentialType << std::endl;
    oss << tabS << "Bond potential params = ["
        << util::io::printStr<double>(d_bondPotentialParams, 0) << "]"
        << std::endl;
    oss << tabS << "State potential type = " << d_statePotentialType
        << std::endl;
    oss << tabS << "State potential params = ["
        << util::io::printStr<double>(d_statePotentialParams, 0) << "]"
        << std::endl;
    oss << tabS << "Influence function type = " << d_influenceFnType
        << std::endl;
    oss << tabS << "Influence function params = ["
        << util::io::printStr<double>(d_influenceFnParams, 0) << "]"
        << std::endl;
    oss << tabS
        << "Irreversible bond breaking enabled = " << d_irreversibleBondBreak
        << std::endl;
    oss << tabS << "State contribution from broken bond enabled = "
        << d_stateContributionFromBrokenBond << std::endl;
    oss << tabS << "Check Sc factor = " << d_checkScFactor << std::endl;
    oss << tabS << "Compute parameters from elastic properties = "
        << d_computeParamsFromElastic << std::endl;
    oss << d_matData.printStr(nt + 1, lvl);
    oss << tabS << "Density = " << d_density << std::endl;
    oss << tabS << "Horizon = " << d_horizon << std::endl;
    oss << tabS << "Horizon to mesh ratio = " << d_horizonMeshRatio
        << std::endl;
    oss << tabS << std::endl;
 
    return oss.str();
  }
 
  void print(int nt = 0, int lvl = 0) const { std::cout << printStr(nt, lvl); }
};
 
} // namespace inp
 
#endif // INP_MATERIALDECK_H