diff --git a/src/GRANULAR/fix_granular_mdr.cpp b/src/GRANULAR/fix_granular_mdr.cpp index b9c245416f..68d9dd4e5e 100644 --- a/src/GRANULAR/fix_granular_mdr.cpp +++ b/src/GRANULAR/fix_granular_mdr.cpp @@ -54,7 +54,15 @@ enum {MEAN_SURF_DISP, RADIUS_UPDATE}; FixGranularMDR::FixGranularMDR(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { - comm_forward = 20; // value needs to match number of values you communicate + comm_forward = 20; // value needs to match number of values you communicate + id_fix = nullptr; +} + +/* ---------------------------------------------------------------------- */ + +FixGranularMDR::~FixGranularMDR() +{ + if (id_fix) modify->delete_fix(id_fix); } /* ---------------------------------------------------------------------- */ @@ -68,13 +76,16 @@ int FixGranularMDR::setmask() /* ---------------------------------------------------------------------- */ -void FixGranularMDR::setup_pre_force(int /*vflag*/) +void FixGranularMDR::post_constructor() { int tmp1, tmp2; + id_fix = utils::strdup("MDR_PARTICLE_HISTORY_VARIABLES"); + modify->add_fix(fmt::format("{} all property/atom d_Ro d_Vcaps d_Vgeo d_Velas d_eps_bar d_dRnumerator d_dRdenominator d_Acon0 d_Acon1 d_Atot d_Atot_sum d_ddelta_bar d_psi d_psi_b d_history_setup_flag d_sigmaxx d_sigmayy d_sigmazz d_contacts d_adhesive_length ghost yes", id_fix)); + int index_Ro = atom->find_custom("Ro",tmp1,tmp2); + int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2); int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2); int index_Velas = atom->find_custom("Velas",tmp1,tmp2); - int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2); int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2); int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2); int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2); @@ -85,12 +96,13 @@ void FixGranularMDR::setup_pre_force(int /*vflag*/) int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); int index_psi = atom->find_custom("psi",tmp1,tmp2); int index_psi_b = atom->find_custom("psi_b",tmp1,tmp2); + int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2); int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2); int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2); int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2); - int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2); int index_contacts = atom->find_custom("contacts",tmp1,tmp2); int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); + Ro = atom->dvector[index_Ro]; Vgeo = atom->dvector[index_Vgeo]; Velas = atom->dvector[index_Velas]; @@ -111,7 +123,12 @@ void FixGranularMDR::setup_pre_force(int /*vflag*/) history_setup_flag = atom->dvector[index_history_setup_flag]; contacts = atom->dvector[index_contacts]; adhesive_length = atom->dvector[index_adhesive_length]; +} +/* ---------------------------------------------------------------------- */ + +void FixGranularMDR::setup_pre_force(int /*vflag*/) +{ pre_force(0); } @@ -119,48 +136,6 @@ void FixGranularMDR::setup_pre_force(int /*vflag*/) void FixGranularMDR::setup(int /*vflag*/) { - int tmp1, tmp2; - int index_Ro = atom->find_custom("Ro",tmp1,tmp2); - int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2); - int index_Velas = atom->find_custom("Velas",tmp1,tmp2); - int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2); - int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2); - int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2); - int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2); - int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); - int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2); - int index_Atot = atom->find_custom("Atot",tmp1,tmp2); - int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2); - int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); - int index_psi = atom->find_custom("psi",tmp1,tmp2); - int index_psi_b = atom->find_custom("psi_b",tmp1,tmp2); - int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2); - int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2); - int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2); - int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2); - int index_contacts = atom->find_custom("contacts",tmp1,tmp2); - int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); - Ro = atom->dvector[index_Ro]; - Vgeo = atom->dvector[index_Vgeo]; - Velas = atom->dvector[index_Velas]; - Vcaps = atom->dvector[index_Vcaps]; - eps_bar = atom->dvector[index_eps_bar]; - dRnumerator = atom->dvector[index_dRnumerator]; - dRdenominator = atom->dvector[index_dRdenominator]; - Acon0 = atom->dvector[index_Acon0]; - Acon1 = atom->dvector[index_Acon1]; - Atot = atom->dvector[index_Atot]; - Atot_sum = atom->dvector[index_Atot_sum]; - ddelta_bar = atom->dvector[index_ddelta_bar]; - psi = atom->dvector[index_psi]; - psi_b = atom->dvector[index_psi_b]; - sigmaxx = atom->dvector[index_sigmaxx]; - sigmayy = atom->dvector[index_sigmayy]; - sigmazz = atom->dvector[index_sigmazz]; - history_setup_flag = atom->dvector[index_history_setup_flag]; - contacts = atom->dvector[index_contacts]; - adhesive_length = atom->dvector[index_adhesive_length]; - end_of_step(); } @@ -381,121 +356,121 @@ void FixGranularMDR::calculate_contact_penalty() double radi = radius[i]; jlist = firstneigh[i]; jnum = numneigh[i]; - for (int jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - jtype = type[j]; - double radj = radius[j]; - const double delx_ij = x[j][0] - xtmp; - const double dely_ij = x[j][1] - ytmp; - const double delz_ij = x[j][2] - ztmp; - const double rsq_ij = delx_ij*delx_ij + dely_ij*dely_ij + delz_ij*delz_ij; - const double r_ij = sqrt(rsq_ij); - const double rinv_ij = 1.0/r_ij; - const double radsum_ij = radi + radj; - const double deltan_ij = radsum_ij - r_ij; - if (deltan_ij >= 0.0) { - for (int kk = jj; kk < jnum; kk++) { - k = jlist[kk]; - k &= NEIGHMASK; - ktype = type[k]; - if (kk != jj) { - const double delx_ik = x[k][0] - xtmp; - const double dely_ik = x[k][1] - ytmp; - const double delz_ik = x[k][2] - ztmp; - const double rsq_ik = delx_ik*delx_ik + dely_ik*dely_ik + delz_ik*delz_ik; - const double r_ik = sqrt(rsq_ik); - const double rinv_ik = 1.0/r_ik; - const double radk = radius[k]; - const double radsum_ik = radi + radk; - const double deltan_ik = radsum_ik - r_ik; - const double delx_jk = x[k][0] - x[j][0]; - const double dely_jk = x[k][1] - x[j][1]; - const double delz_jk = x[k][2] - x[j][2]; - const double rsq_jk = delx_jk*delx_jk + dely_jk*dely_jk + delz_jk*delz_jk; - const double r_jk = sqrt(rsq_jk); - const double rinv_jk = 1.0/r_jk; - const double radsum_jk = radj + radk; - const double deltan_jk = radsum_jk - r_jk; - if (deltan_ik >= 0.0 && deltan_jk >= 0.0) { + for (int jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtype = type[j]; + double radj = radius[j]; + const double delx_ij = x[j][0] - xtmp; + const double dely_ij = x[j][1] - ytmp; + const double delz_ij = x[j][2] - ztmp; + const double rsq_ij = delx_ij*delx_ij + dely_ij*dely_ij + delz_ij*delz_ij; + const double r_ij = sqrt(rsq_ij); + const double rinv_ij = 1.0/r_ij; + const double radsum_ij = radi + radj; + const double deltan_ij = radsum_ij - r_ij; + if (deltan_ij >= 0.0) { + for (int kk = jj; kk < jnum; kk++) { + k = jlist[kk]; + k &= NEIGHMASK; + ktype = type[k]; + if (kk != jj) { + const double delx_ik = x[k][0] - xtmp; + const double dely_ik = x[k][1] - ytmp; + const double delz_ik = x[k][2] - ztmp; + const double rsq_ik = delx_ik*delx_ik + dely_ik*dely_ik + delz_ik*delz_ik; + const double r_ik = sqrt(rsq_ik); + const double rinv_ik = 1.0/r_ik; + const double radk = radius[k]; + const double radsum_ik = radi + radk; + const double deltan_ik = radsum_ik - r_ik; + const double delx_jk = x[k][0] - x[j][0]; + const double dely_jk = x[k][1] - x[j][1]; + const double delz_jk = x[k][2] - x[j][2]; + const double rsq_jk = delx_jk*delx_jk + dely_jk*dely_jk + delz_jk*delz_jk; + const double r_jk = sqrt(rsq_jk); + const double rinv_jk = 1.0/r_jk; + const double radsum_jk = radj + radk; + const double deltan_jk = radsum_jk - r_jk; + if (deltan_ik >= 0.0 && deltan_jk >= 0.0) { - // pull ij history - history_ij = &allhistory[size_history * jj]; - double * pij = &history_ij[22]; // penalty for contact i and j + // pull ij history + history_ij = &allhistory[size_history * jj]; + double * pij = &history_ij[22]; // penalty for contact i and j - // pull ik history - history_ik = &allhistory[size_history * kk]; - double * pik = &history_ik[22]; // penalty for contact i and k + // pull ik history + history_ik = &allhistory[size_history * kk]; + double * pik = &history_ik[22]; // penalty for contact i and k - // we don't know if who owns the contact ahead of time, k might be in j's neigbor list or vice versa, - // so we need to manually search to figure out the owner check if k is in the neighbor list of j - double * pjk = NULL; - int * const jklist = firstneigh[j]; - const int jknum = numneigh[j]; - for (int jk = 0; jk < jknum; jk++) { - const int kneigh = jklist[jk] & NEIGHMASK; - if (k == kneigh) { - allhistory_j = firsthistory[j]; - history_jk = &allhistory_j[size_history * jk]; - pjk = &history_jk[22]; // penalty for contact j and k + // we don't know if who owns the contact ahead of time, k might be in j's neigbor list or vice versa, + // so we need to manually search to figure out the owner check if k is in the neighbor list of j + double * pjk = NULL; + int * const jklist = firstneigh[j]; + const int jknum = numneigh[j]; + for (int jk = 0; jk < jknum; jk++) { + const int kneigh = jklist[jk] & NEIGHMASK; + if (k == kneigh) { + allhistory_j = firsthistory[j]; + history_jk = &allhistory_j[size_history * jk]; + pjk = &history_jk[22]; // penalty for contact j and k + break; + } + } + + // check if j is in the neighbor list of k + if (pjk == NULL) { + int * const kjlist = firstneigh[k]; + const int kjnum = numneigh[k]; + for (int kj = 0; kj < kjnum; kj++) { + const int jneigh = kjlist[kj] & NEIGHMASK; + if (j == jneigh) { + allhistory_k = firsthistory[k]; + history_kj = &allhistory_k[size_history * kj]; + pjk = &history_kj[22]; // penalty for contact j and k break; } } + } - // check if j is in the neighbor list of k - if (pjk == NULL) { - int * const kjlist = firstneigh[k]; - const int kjnum = numneigh[k]; - for (int kj = 0; kj < kjnum; kj++) { - const int jneigh = kjlist[kj] & NEIGHMASK; - if (j == jneigh) { - allhistory_k = firsthistory[k]; - history_kj = &allhistory_k[size_history * kj]; - pjk = &history_kj[22]; // penalty for contact j and k - break; - } - } - } - - std::vector distances = {r_ij,r_ik,r_jk}; - auto maxElement = std::max_element(distances.begin(), distances.end()); - double maxValue = *maxElement; - int maxIndex = std::distance(distances.begin(), maxElement); - if (maxIndex == 0) { // the central particle is k - const double enx_ki = -delx_ik * rinv_ik; - const double eny_ki = -dely_ik * rinv_ik; - const double enz_ki = -delz_ik * rinv_ik; - const double enx_kj = -delx_jk * rinv_jk; - const double eny_kj = -dely_jk * rinv_jk; - const double enz_kj = -delz_jk * rinv_jk; - const double alpha = std::acos(enx_ki*enx_kj + eny_ki*eny_kj + enz_ki*enz_kj); - pij[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0)) ); - } else if (maxIndex == 1) { // the central particle is j - const double enx_ji = -delx_ij * rinv_ij; - const double eny_ji = -dely_ij * rinv_ij; - const double enz_ji = -delz_ij * rinv_ij; - const double enx_jk = delx_jk * rinv_jk; - const double eny_jk = dely_jk * rinv_jk; - const double enz_jk = delz_jk * rinv_jk; - const double alpha = std::acos(enx_ji*enx_jk + eny_ji*eny_jk + enz_ji*enz_jk); - pik[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0)) ); - } else { // the central particle is i - if (j < atom->nlocal || k < atom->nlocal) { - const double enx_ij = delx_ij * rinv_ij; - const double eny_ij = dely_ij * rinv_ij; - const double enz_ij = delz_ij * rinv_ij; - const double enx_ik = delx_ik * rinv_ik; - const double eny_ik = dely_ik * rinv_ik; - const double enz_ik = delz_ik * rinv_ik; - const double alpha = std::acos(enx_ij*enx_ik + eny_ij*eny_ik + enz_ij*enz_ik); - pjk[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0)) ); - } + std::vector distances = {r_ij,r_ik,r_jk}; + auto maxElement = std::max_element(distances.begin(), distances.end()); + double maxValue = *maxElement; + int maxIndex = std::distance(distances.begin(), maxElement); + if (maxIndex == 0) { // the central particle is k + const double enx_ki = -delx_ik * rinv_ik; + const double eny_ki = -dely_ik * rinv_ik; + const double enz_ki = -delz_ik * rinv_ik; + const double enx_kj = -delx_jk * rinv_jk; + const double eny_kj = -dely_jk * rinv_jk; + const double enz_kj = -delz_jk * rinv_jk; + const double alpha = std::acos(enx_ki*enx_kj + eny_ki*eny_kj + enz_ki*enz_kj); + pij[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0))); + } else if (maxIndex == 1) { // the central particle is j + const double enx_ji = -delx_ij * rinv_ij; + const double eny_ji = -dely_ij * rinv_ij; + const double enz_ji = -delz_ij * rinv_ij; + const double enx_jk = delx_jk * rinv_jk; + const double eny_jk = dely_jk * rinv_jk; + const double enz_jk = delz_jk * rinv_jk; + const double alpha = std::acos(enx_ji*enx_jk + eny_ji*eny_jk + enz_ji*enz_jk); + pik[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0))); + } else { // the central particle is i + if (j < atom->nlocal || k < atom->nlocal) { + const double enx_ij = delx_ij * rinv_ij; + const double eny_ij = dely_ij * rinv_ij; + const double enz_ij = delz_ij * rinv_ij; + const double enx_ik = delx_ik * rinv_ik; + const double eny_ik = dely_ik * rinv_ik; + const double enz_ik = delz_ik * rinv_ik; + const double alpha = std::acos(enx_ij*enx_ik + eny_ij*eny_ik + enz_ij*enz_ik); + pjk[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0))); } } } } } } + } } } @@ -511,16 +486,12 @@ void FixGranularMDR::mean_surf_disp() NeighList * list = pair->list; const int size_history = pair->get_size_history(); - - int i,j,k,ii,jj,inum,jnum,itype,jtype; - int *ilist,*jlist,*numneigh,**firstneigh; int *touch,**firsttouch; double *history,*allhistory,**firsthistory; bool touchflag = false; - class GranularModel* model; class GranularModel** models_list = pair->models_list; int ** types_indices = pair->types_indices; @@ -530,7 +501,6 @@ void FixGranularMDR::mean_surf_disp() double *radius = atom->radius; int nlocal = atom->nlocal; - inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -657,11 +627,9 @@ void FixGranularMDR::update_fix_gran_wall(Fix* fix_in) Region * region = fix->region; const int size_history = model->size_history; - int i, m, nc, iwall; double vwall[3]; bool touchflag = false; - int history_update = 1; model->history_update = history_update; @@ -688,21 +656,21 @@ void FixGranularMDR::update_fix_gran_wall(Fix* fix_in) nc = region->surface(x[i][0], x[i][1], x[i][2], radius[i] + model->pulloff_distance(radius[i], 0.0)); - if (nc == 0) { - fix->ncontact[i] = 0; - continue; - } - if (nc == 1) { - fix->c2r[0] = 0; - iwall = region->contact[0].iwall; - if (fix->ncontact[i] == 0) { - fix->ncontact[i] = 1; - fix->walls[i][0] = iwall; - for (m = 0; m < size_history; m++) fix->history_many[i][0][m] = 0.0; - } else if (fix->ncontact[i] > 1 || iwall != fix->walls[i][0]) - fix->update_contacts(i, nc); - } else + if (nc == 0) { + fix->ncontact[i] = 0; + continue; + } + if (nc == 1) { + fix->c2r[0] = 0; + iwall = region->contact[0].iwall; + if (fix->ncontact[i] == 0) { + fix->ncontact[i] = 1; + fix->walls[i][0] = iwall; + for (m = 0; m < size_history; m++) fix->history_many[i][0][m] = 0.0; + } else if (fix->ncontact[i] > 1 || iwall != fix->walls[i][0]) fix->update_contacts(i, nc); + } else + fix->update_contacts(i, nc); // process current contacts diff --git a/src/GRANULAR/fix_granular_mdr.h b/src/GRANULAR/fix_granular_mdr.h index 149bb61ec7..319300cc7b 100644 --- a/src/GRANULAR/fix_granular_mdr.h +++ b/src/GRANULAR/fix_granular_mdr.h @@ -26,44 +26,46 @@ namespace LAMMPS_NS { class FixGranularMDR : public Fix { public: - double * Ro; - double * Vgeo; - double * Velas; - double * Vcaps; - double * eps_bar; - double * dRnumerator; - double * dRdenominator; - double * Acon0; - double * Acon1; - double * Atot; - double * Atot_sum; - double * ddelta_bar; - double * psi; - double * psi_b; - double * sigmaxx; - double * sigmayy; - double * sigmazz; - double * history_setup_flag; - double * contacts; - double * adhesive_length; - FixGranularMDR(class LAMMPS *, int, char **); + ~FixGranularMDR() override; int setmask() override; + void post_constructor() override; void setup(int) override; void setup_pre_force(int) override; void pre_force(int) override; - void end_of_step() override; // FOR MDR + void end_of_step() override; int pack_forward_comm(int, int *, double *, int, int *) override; void unpack_forward_comm(int, int, double *) override; private: int comm_stage; + char *id_fix; void radius_update(); void mean_surf_disp(); void calculate_contact_penalty(); void update_fix_gran_wall(Fix*); + double *Ro; // initial radius + double *Vgeo; // geometric particle volume of apparent particle afterremoving spherical cap volume + double *Velas; // particle volume from linear elasticity + double *Vcaps; // spherical cap volume from intersection of apparentradius particle and contact planes + double *eps_bar; // volume-averaged infinitesimal strain tensor + double *dRnumerator; // summation of numerator terms in calculation of dR + double *dRdenominator; // summation of denominator terms in calculation of dR + double *Acon0; // total area involved in contacts: Acon^{n} + double *Acon1; // total area involved in contacts: Acon^{n+1} + double *Atot; // total particle area + double *Atot_sum; // running sum of contact area minus cap area + double *ddelta_bar; // change in mean surface displacement + double *psi; // ratio of free surface area to total surface area + double *psi_b; // TEMPORARY, SINCE PSI_B IS ALREADY DEFINED IN THEINPUT SCRIPT + double *sigmaxx; // xx-component of the stress tensor, not necessary forforce calculation + double *sigmayy; // yy-component of the stress tensor, not necessary forforce calculation + double *sigmazz; // zz-component of the stress tensor, not necessary forforce calculation + double *history_setup_flag; // flag to check if history variables have beeninitialized + double *contacts; // total contacts on particle + double *adhesive_length; // total length of adhesive contact on a particle }; } // namespace LAMMPS_NS diff --git a/src/GRANULAR/gran_sub_mod_normal.cpp b/src/GRANULAR/gran_sub_mod_normal.cpp index 43e77b1c4e..9eddc827ff 100644 --- a/src/GRANULAR/gran_sub_mod_normal.cpp +++ b/src/GRANULAR/gran_sub_mod_normal.cpp @@ -13,12 +13,13 @@ #include "gran_sub_mod_normal.h" +#include "atom.h" #include "error.h" +#include "citeme.h" #include "granular_model.h" #include "math_const.h" -#include "atom.h" +#include "modify.h" #include "update.h" -#include "citeme.h" #include #include @@ -414,6 +415,7 @@ GranSubModNormalMDR::GranSubModNormalMDR(GranularModel *gm, LAMMPS *lmp) : num_coeffs = 6; // Young's Modulus, Poisson's ratio, yield stress, effective surface energy, psi_b, coefficent of restitution contact_radius_flag = 1; size_history = 26; + fix_mdr_flag = 0; nondefault_history_transfer = 1; transfer_history_factor = new double[size_history]; @@ -425,6 +427,14 @@ GranSubModNormalMDR::GranSubModNormalMDR(GranularModel *gm, LAMMPS *lmp) : /* ---------------------------------------------------------------------- */ +GranSubModNormalMDR::~GranSubModNormalMDR() +{ + if (fix_mdr_flag) + modify->delete_fix("MDR"); +} + +/* ---------------------------------------------------------------------- */ + void GranSubModNormalMDR::coeffs_to_local() { E = coeffs[0]; // Young's modulus @@ -444,6 +454,57 @@ void GranSubModNormalMDR::coeffs_to_local() /* ---------------------------------------------------------------------- */ +void GranSubModNormalMDR::init() +{ + if (!fix_mdr_flag) { + if (modify->get_fix_by_style("MDR").size() == 0) + modify->add_fix("MDR all GRANULAR/MDR"); + fix_mdr_flag = 1; + } + + // initialize particle history variables + int tmp1, tmp2; + int index_Ro = atom->find_custom("Ro", tmp1, tmp2); // initial radius + int index_Vcaps = atom->find_custom("Vcaps", tmp1, tmp2); // spherical cap volume from intersection of apparent radius particle and contact planes + int index_Vgeo = atom->find_custom("Vgeo", tmp1, tmp2); // geometric particle volume of apparent particle after removing spherical cap volume + int index_Velas = atom->find_custom("Velas", tmp1, tmp2); // particle volume from linear elasticity + int index_eps_bar = atom->find_custom("eps_bar", tmp1, tmp2); // volume-averaged infinitesimal strain tensor + int index_dRnumerator = atom->find_custom("dRnumerator", tmp1, tmp2); // summation of numerator terms in calculation of dR + int index_dRdenominator = atom->find_custom("dRdenominator", tmp1, tmp2); // summation of denominator terms in calculation of dR + int index_Acon0 = atom->find_custom("Acon0", tmp1, tmp2); // total area involved in contacts: Acon^{n} + int index_Acon1 = atom->find_custom("Acon1", tmp1, tmp2); // total area involved in contacts: Acon^{n+1} + int index_Atot = atom->find_custom("Atot", tmp1, tmp2); // total particle area + int index_Atot_sum = atom->find_custom("Atot_sum", tmp1, tmp2); // running sum of contact area minus cap area + int index_ddelta_bar = atom->find_custom("ddelta_bar", tmp1, tmp2); // change in mean surface displacement + int index_psi = atom->find_custom("psi", tmp1, tmp2); // ratio of free surface area to total surface area + int index_sigmaxx = atom->find_custom("sigmaxx", tmp1, tmp2); // xx-component of the stress tensor, not necessary for force calculation + int index_sigmayy = atom->find_custom("sigmayy", tmp1, tmp2); // yy-component of the stress tensor, not necessary for force calculation + int index_sigmazz = atom->find_custom("sigmazz", tmp1, tmp2); // zz-component of the stress tensor, not necessary for force calculation + int index_contacts = atom->find_custom("contacts", tmp1, tmp2); // total contacts on particle + int index_adhesive_length = atom->find_custom("adhesive_length", tmp1, tmp2); // total contacts on particle + + Rinitial = atom->dvector[index_Ro]; + Vgeo = atom->dvector[index_Vgeo]; + Velas = atom->dvector[index_Velas]; + Vcaps = atom->dvector[index_Vcaps]; + eps_bar = atom->dvector[index_eps_bar]; + dRnumerator = atom->dvector[index_dRnumerator]; + dRdenominator = atom->dvector[index_dRdenominator]; + Acon0 = atom->dvector[index_Acon0]; + Acon1 = atom->dvector[index_Acon1]; + Atot = atom->dvector[index_Atot]; + Atot_sum = atom->dvector[index_Atot_sum]; + ddelta_bar = atom->dvector[index_ddelta_bar]; + psi = atom->dvector[index_psi]; + sigmaxx = atom->dvector[index_sigmaxx]; + sigmayy = atom->dvector[index_sigmayy]; + sigmazz = atom->dvector[index_sigmazz]; + contacts = atom->dvector[index_contacts]; + adhesive_length = atom->dvector[index_adhesive_length]; +} + +/* ---------------------------------------------------------------------- */ + double GranSubModNormalMDR::calculate_forces() { @@ -507,46 +568,6 @@ double GranSubModNormalMDR::calculate_forces() const int deltap_offset_0 = 24; const int deltap_offset_1 = 25; - - // initialize particle history variables - int tmp1, tmp2; - int index_Ro = atom->find_custom("Ro",tmp1,tmp2); // initial radius - int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2); // spherical cap volume from intersection of apparent radius particle and contact planes - int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2); // geometric particle volume of apparent particle after removing spherical cap volume - int index_Velas = atom->find_custom("Velas",tmp1,tmp2); // particle volume from linear elasticity - int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2); // volume-averaged infinitesimal strain tensor - int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2); // summation of numerator terms in calculation of dR - int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2); // summation of denominator terms in calculation of dR - int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); // total area involved in contacts: Acon^{n} - int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2); // total area involved in contacts: Acon^{n+1} - int index_Atot = atom->find_custom("Atot",tmp1,tmp2); // total particle area - int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2); // running sum of contact area minus cap area - int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); // change in mean surface displacement - int index_psi = atom->find_custom("psi",tmp1,tmp2); // ratio of free surface area to total surface area - int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2); // xx-component of the stress tensor, not necessary for force calculation - int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2); // yy-component of the stress tensor, not necessary for force calculation - int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2); // zz-component of the stress tensor, not necessary for force calculation - int index_contacts = atom->find_custom("contacts",tmp1,tmp2); // total contacts on particle - int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); // total contacts on particle - double * Rinitial = atom->dvector[index_Ro]; - double * Vgeo = atom->dvector[index_Vgeo]; - double * Velas = atom->dvector[index_Velas]; - double * Vcaps = atom->dvector[index_Vcaps]; - double * eps_bar = atom->dvector[index_eps_bar]; - double * dRnumerator = atom->dvector[index_dRnumerator]; - double * dRdenominator = atom->dvector[index_dRdenominator]; - double * Acon0 = atom->dvector[index_Acon0]; - double * Acon1 = atom->dvector[index_Acon1]; - double * Atot = atom->dvector[index_Atot]; - double * Atot_sum = atom->dvector[index_Atot_sum]; - double * ddelta_bar = atom->dvector[index_ddelta_bar]; - double * psi = atom->dvector[index_psi]; - double * sigmaxx = atom->dvector[index_sigmaxx]; - double * sigmayy = atom->dvector[index_sigmayy]; - double * sigmazz = atom->dvector[index_sigmazz]; - double * contacts = atom->dvector[index_contacts]; - double * adhesive_length = atom->dvector[index_adhesive_length]; - double * history = & gm->history[history_index]; // load in all history variables // Rigid flat placement scheme diff --git a/src/GRANULAR/gran_sub_mod_normal.h b/src/GRANULAR/gran_sub_mod_normal.h index bfd6121ec8..29e201245a 100644 --- a/src/GRANULAR/gran_sub_mod_normal.h +++ b/src/GRANULAR/gran_sub_mod_normal.h @@ -139,13 +139,20 @@ namespace Granular_NS { class GranSubModNormalMDR : public GranSubModNormal { public: GranSubModNormalMDR(class GranularModel *, class LAMMPS *); + ~GranSubModNormalMDR() override; void coeffs_to_local() override; + void init() override; double calculate_forces() override; void set_fncrit() override; double psi_b; protected: double E, nu, Y, gamma, CoR, F; + + double *Rinitial, *Vgeo, *Velas, *Vcaps, *eps_bar, *dRnumerator; + double *dRdenominator, *Acon0, *Acon1, *Atot, *Atot_sum, *ddelta_bar; + double *psi, *sigmaxx, *sigmayy, *sigmazz, *contacts, *adhesive_length; + int fix_mdr_flag; }; } // namespace Granular_NS diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 23522608c8..fe0c75a2af 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -78,8 +78,6 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) fix_history = nullptr; fix_dummy = dynamic_cast(modify->add_fix("NEIGH_HISTORY_GRANULAR_DUMMY all DUMMY")); - - fix_flag = 0; } /* ---------------------------------------------------------------------- */ @@ -475,41 +473,6 @@ void PairGranular::init_style() if (!fix_history) error->all(FLERR,"Could not find pair fix neigh history ID"); } - if (model->sub_models[NORMAL]->name == "mdr") { - //Store persistent per atom quantities - if (! fix_flag) { - int tmp1, tmp2; - const char * id_fix = "MDR_PARTICLE_HISTORY_VARIABLES"; - modify->add_fix(fmt::format("{} all property/atom d_Ro d_Vcaps d_Vgeo d_Velas d_eps_bar d_dRnumerator d_dRdenominator d_Acon0 d_Acon1 d_Atot d_Atot_sum d_ddelta_bar d_psi d_psi_b d_history_setup_flag d_sigmaxx d_sigmayy d_sigmazz d_contacts d_adhesive_length ghost yes", id_fix)); - // d2_volSums 4 --> allows an array of 4 to defined. - index_Ro = atom->find_custom("Ro",tmp1,tmp2); // initial radius - index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2); // spherical cap volume from intersection of apparent radius particle and contact planes - index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2); // geometric particle volume of apparent particle after removing spherical cap volume - index_Velas = atom->find_custom("Velas",tmp1,tmp2); // particle volume from linear elasticity - index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2); // volume-averaged infinitesimal strain tensor - index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2); // summation of numerator terms in calculation of dR - index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2); // summation of denominator terms in calculation of dR - index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); // total area involved in contacts: Acon^{n} - index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2); // total area involved in contacts: Acon^{n+1} - index_Atot = atom->find_custom("Atot",tmp1,tmp2); // total particle area - index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2); // running sum of contact area minus cap area - index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); // change in mean surface displacement - index_psi = atom->find_custom("psi",tmp1,tmp2); // ratio of free surface area to total surface area - index_psi_b = atom->find_custom("psi_b",tmp1,tmp2); // TEMPORARY, SINCE PSI_B IS ALREADY DEFINED IN THE INPUT SCRIPT - index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2); // flag to check if history variables have been initialized - index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2); // xx-component of the stress tensor, not necessary for force calculation - index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2); // yy-component of the stress tensor, not necessary for force calculation - index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2); // zz-component of the stress tensor, not necessary for force calculation - index_contacts = atom->find_custom("contacts",tmp1,tmp2); // total contacts on particle - index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); // total length of adhesive contact on a particle - - // Initiate MDR radius update fix - modify->add_fix("fix_mdr all GRANULAR/MDR"); - - fix_flag = 1; - } - } - // check for FixFreeze and set freeze_group_bit auto fixlist = modify->get_fix_by_style("^freeze"); diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 490133be09..1454291cd8 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -53,8 +53,6 @@ class PairGranular : public Pair { int **types_indices; int nmodels, maxmodels; - class FixStoreAtom * fix_store; - protected: int freeze_group_bit; int use_history; @@ -67,29 +65,6 @@ class PairGranular : public Pair { class FixDummy *fix_dummy; class FixNeighHistory *fix_history; - // MDR particle history variables - int fix_flag; - int index_Ro; - int index_Vcaps; - int index_Vgeo; - int index_Velas; - int index_eps_bar; - int index_dRnumerator; - int index_dRdenominator; - int index_Acon0; - int index_Acon1; - int index_Atot; - int index_Atot_sum; - int index_ddelta_bar; - int index_psi; - int index_psi_b; - int index_history_setup_flag; - int index_sigmaxx; - int index_sigmayy; - int index_sigmazz; - int index_contacts; - int index_adhesive_length; - // storage of rigid body masses for use in granular interactions class Fix *fix_rigid; // ptr to rigid body fix, null pointer if none