#include "gran_sub_mod_custom.h"
#include "gran_sub_mod_normal.h"
#include "granular_model.h"
#include "error.h" // 确保包含错误头文件

using namespace LAMMPS_NS;
using namespace Granular_NS;

static constexpr double FOURTHIRDS = 1.33333333333333333333; // 4/3

GranSubModNormalHertzPiecewise::GranSubModNormalHertzPiecewise(GranularModel *gm, LAMMPS *lmp)
    : GranSubModNormal(gm, lmp) {
  material_properties = 1;
  num_coeffs = 6;
  contact_radius_flag = 1;
  mixed_coefficients = 0;  // 初始化 mixed_coefficients
}

void GranSubModNormalHertzPiecewise::coeffs_to_local() {
  Emod1 = coeffs[0];
  Emod2 = coeffs[1];
  poiss1 = coeffs[2];
  poiss2 = coeffs[3];
  damp = coeffs[4];
  delta_switch = coeffs[5];

  if (Emod1 < 0.0 || Emod2 < 0.0 || damp < 0.0) {
    error->all(FLERR, "Illegal Hertz material normal model");
  }

}


double GranSubModNormalHertzPiecewise::calculate_forces() {
  double Fne;

  // 当接触位移大于阈值时，使用不同的刚度计算力
  if (gm->delta >= delta_switch) {
    Fne = gm->contact_radius * (FOURTHIRDS * mix_stiffnessE(Emod1, Emod1, poiss1, poiss1) * delta_switch +
                                FOURTHIRDS * mix_stiffnessE(Emod2, Emod2, poiss2, poiss2) * (gm->delta - delta_switch));
  } else {
    Fne = FOURTHIRDS * mix_stiffnessE(Emod1, Emod1, poiss1, poiss1) * gm->contact_radius * gm->delta;  // 小于阈值时使用k1
  }

  return Fne;
}