From f51ab2c4401c9fe6e01e7b3cb72bf97bacb4b200 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 16 Nov 2022 16:36:22 -0700 Subject: [PATCH] Restarting limit_damping, other minor simplifications/cleanups --- examples/granular/in.pour.drum | 2 +- src/GRANULAR/fix_wall_gran.cpp | 32 ++++++++--------- src/GRANULAR/fix_wall_gran.h | 3 +- src/GRANULAR/fix_wall_gran_region.cpp | 2 +- src/GRANULAR/gran_sub_mod.h | 4 +-- src/GRANULAR/gran_sub_mod_damping.cpp | 2 +- src/GRANULAR/gran_sub_mod_damping.h | 2 -- src/GRANULAR/gran_sub_mod_heat.h | 2 -- src/GRANULAR/gran_sub_mod_normal.h | 4 +-- src/GRANULAR/gran_sub_mod_rolling.h | 2 -- src/GRANULAR/gran_sub_mod_tangential.h | 4 --- src/GRANULAR/gran_sub_mod_twisting.h | 4 +-- src/GRANULAR/granular_model.cpp | 44 ++++++++++++----------- src/GRANULAR/granular_model.h | 8 ++--- src/GRANULAR/pair_granular.cpp | 50 +++++++++++--------------- 15 files changed, 72 insertions(+), 93 deletions(-) diff --git a/examples/granular/in.pour.drum b/examples/granular/in.pour.drum index e0a0455f61..6c93b700d6 100644 --- a/examples/granular/in.pour.drum +++ b/examples/granular/in.pour.drum @@ -44,7 +44,7 @@ change_box all boundary p p f pair_style granular pair_coeff 1 * hertz/material 1e5 0.2 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji -pair_coeff 2 2 jkr 1e5 0.1 0.3 50 tangential mindlin NULL 1.0 0.5 rolling sds 1e3 1e3 0.1 twisting marshall damping tsuji +pair_coeff 2 2 jkr 1e5 0.1 0.3 50 tangential mindlin NULL 1.0 0.5 rolling sds 1e3 1e3 0.1 twisting marshall variable theta equal 0 diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index fe2286a0fa..3afb0dc60e 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -67,11 +67,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : model = new GranularModel(lmp); model->contact_type = WALL; - if (strcmp(arg[3],"granular") == 0) classic_flag = 0; - else classic_flag = 1; - limit_damping = 0; heat_flag = 0; - int Twall_defined = 0; + int classic_flag = 1; + if (strcmp(arg[3],"granular") == 0) classic_flag = 0; // wall/particle coefficients @@ -117,10 +115,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : } } - // Define default damping model if unspecified, takes no args + // Define default damping submodel if unspecified, takes no args if (!model->damping_model) model->construct_submodel("viscoelastic", DAMPING); - model->init(); size_history = model->size_history; @@ -177,6 +174,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : wiggle = 0; wshear = 0; peratom_flag = 0; + int Twall_defined = 0; while (iarg < narg) { if (strcmp(arg[iarg],"wiggle") == 0) { @@ -216,7 +214,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : } if (heat_flag != Twall_defined) - error->all(FLERR, "To model conduction, must define both heat model and wall temperature"); + error->all(FLERR, "Must define wall temperature with heat model"); if (wallstyle == NOSTYLE) error->all(FLERR,"No wall style defined"); @@ -376,6 +374,7 @@ void FixWallGran::post_force(int /*vflag*/) history_update = 1; if (update->setupflag) history_update = 0; + model->history_update = history_update; // if just reneighbored: // update rigid body masses for owned atoms if using FixRigid @@ -407,10 +406,10 @@ void FixWallGran::post_force(int /*vflag*/) if (wiggle) { double arg = omega * (update->ntimestep - time_origin) * dt; if (wallstyle == axis) { - wlo = lo + amplitude - amplitude*cos(arg); - whi = hi + amplitude - amplitude*cos(arg); + wlo = lo + amplitude - amplitude * cos(arg); + whi = hi + amplitude - amplitude * cos(arg); } - vwall[axis] = amplitude*omega*sin(arg); + vwall[axis] = amplitude * omega * sin(arg); } else if (wshear) vwall[axis] = vshear; // loop over all my atoms @@ -473,19 +472,19 @@ void FixWallGran::post_force(int /*vflag*/) if (del1 < del2) dz = del1; else dz = -del2; } else if (wallstyle == ZCYLINDER) { - delxy = sqrt(x[i][0]*x[i][0] + x[i][1]*x[i][1]); + delxy = sqrt(x[i][0] * x[i][0] + x[i][1] * x[i][1]); delr = cylradius - delxy; if (delr > radius[i]) { dz = cylradius; rwall = 0.0; } else { - dx = -delr/delxy * x[i][0]; - dy = -delr/delxy * x[i][1]; + dx = -delr / delxy * x[i][0]; + dy = -delr / delxy * x[i][1]; // rwall = -2r_c if inside cylinder, 2r_c outside - rwall = (delxy < cylradius) ? -2*cylradius : 2*cylradius; + rwall = (delxy < cylradius) ? -2 * cylradius : 2 * cylradius; if (wshear && axis != 2) { - vwall[0] += vshear * x[i][1]/delxy; - vwall[1] += -vshear * x[i][0]/delxy; + vwall[0] += vshear * x[i][1] / delxy; + vwall[1] += -vshear * x[i][0] / delxy; vwall[2] = 0.0; } } @@ -521,7 +520,6 @@ void FixWallGran::post_force(int /*vflag*/) model->meff = meff; model->vi = v[i]; model->omegai = omega[i]; - model->history_update = history_update; if (use_history) model->history = history_one[i]; if (heat_flag) model->Ti = temperature[i]; diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index 16a2fe5dc0..45e4e43844 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -53,7 +53,7 @@ class FixWallGran : public Fix { protected: int wallstyle, wiggle, wshear, axis; - int classic_flag, nlevels_respa; + int nlevels_respa; bigint time_origin; // for granular model choices @@ -69,7 +69,6 @@ class FixWallGran : public Fix { int history_update; // flag for whether shear history is updated int size_history; // # of shear history values per contact int heat_flag; - int limit_damping; int tvar; char *tstr; diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index 7ee794133c..f7b772cb7f 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/src/GRANULAR/fix_wall_gran_region.cpp @@ -135,6 +135,7 @@ void FixWallGranRegion::post_force(int /*vflag*/) history_update = 1; if (update->setupflag) history_update = 0; + model->history_update = history_update; // if just reneighbored: // update rigid body masses for owned atoms if using FixRigid @@ -262,7 +263,6 @@ void FixWallGranRegion::post_force(int /*vflag*/) model->meff = meff; model->vi = v[i]; model->omegai = omega[i]; - model->history_update = history_update; if (use_history) model->history = history_many[i][c2r[ic]]; if (heat_flag) model->Ti = temperature[i]; diff --git a/src/GRANULAR/gran_sub_mod.h b/src/GRANULAR/gran_sub_mod.h index b2335e8cfb..22f1f4e85a 100644 --- a/src/GRANULAR/gran_sub_mod.h +++ b/src/GRANULAR/gran_sub_mod.h @@ -29,8 +29,8 @@ class GranSubMod : protected Pointers { double *coeffs; void read_restart(); virtual void mix_coeffs(double*, double*); - virtual void coeffs_to_local() = 0; - virtual void init() = 0; // called after all submodel coeffs defined + virtual void coeffs_to_local() {}; + virtual void init() {}; // called after all submodels + coeffs defined void allocate_coeffs(); std::string name; diff --git a/src/GRANULAR/gran_sub_mod_damping.cpp b/src/GRANULAR/gran_sub_mod_damping.cpp index e78233735a..820d552679 100644 --- a/src/GRANULAR/gran_sub_mod_damping.cpp +++ b/src/GRANULAR/gran_sub_mod_damping.cpp @@ -108,7 +108,7 @@ void GranSubModDampingTsuji::init() { double tmp = gm->normal_model->damp; damp = 1.2728 - 4.2783 * tmp + 11.087 * square(tmp); - damp += -22.348 * cube(tmp)+ 27.467 * powint(tmp, 4); + damp += -22.348 * cube(tmp) + 27.467 * powint(tmp, 4); damp += -18.022 * powint(tmp, 5) + 4.8218 * powint(tmp,6); } diff --git a/src/GRANULAR/gran_sub_mod_damping.h b/src/GRANULAR/gran_sub_mod_damping.h index e5f4227ce8..5d245c384d 100644 --- a/src/GRANULAR/gran_sub_mod_damping.h +++ b/src/GRANULAR/gran_sub_mod_damping.h @@ -48,8 +48,6 @@ class GranSubModDamping : public GranSubMod { public: GranSubModDamping(class GranularModel *, class LAMMPS *); ~GranSubModDamping() {}; - virtual void coeffs_to_local() {}; - virtual void mix_coeffs(double*, double*) {}; virtual void init(); virtual double calculate_forces() = 0; double damp_prefactor; // Used by tangential models diff --git a/src/GRANULAR/gran_sub_mod_heat.h b/src/GRANULAR/gran_sub_mod_heat.h index a5408a0fc2..ca7e82313a 100644 --- a/src/GRANULAR/gran_sub_mod_heat.h +++ b/src/GRANULAR/gran_sub_mod_heat.h @@ -35,8 +35,6 @@ class GranSubModHeat : public GranSubMod { public: GranSubModHeat(class GranularModel *, class LAMMPS *); ~GranSubModHeat() {}; - virtual void coeffs_to_local() {}; - virtual void init() {}; virtual double calculate_heat() = 0; }; diff --git a/src/GRANULAR/gran_sub_mod_normal.h b/src/GRANULAR/gran_sub_mod_normal.h index 49e7d012c4..7dd9f82349 100644 --- a/src/GRANULAR/gran_sub_mod_normal.h +++ b/src/GRANULAR/gran_sub_mod_normal.h @@ -51,8 +51,6 @@ class GranSubModNormal : public GranSubMod { public: GranSubModNormal(class GranularModel *, class LAMMPS *); ~GranSubModNormal() {}; - virtual void coeffs_to_local() {}; - virtual void init() {}; virtual bool touch(); virtual double pulloff_distance(double, double); virtual double calculate_area(); @@ -60,7 +58,7 @@ class GranSubModNormal : public GranSubMod { virtual void set_fncrit(); double damp; // argument historically needed by damping double Emod, poiss; - double Fncrit, knfac; + double Fncrit; int material_properties, cohesive_flag; }; diff --git a/src/GRANULAR/gran_sub_mod_rolling.h b/src/GRANULAR/gran_sub_mod_rolling.h index 5739067f50..980daeca89 100644 --- a/src/GRANULAR/gran_sub_mod_rolling.h +++ b/src/GRANULAR/gran_sub_mod_rolling.h @@ -35,8 +35,6 @@ class GranSubModRolling : public GranSubMod { public: GranSubModRolling(class GranularModel *, class LAMMPS *); ~GranSubModRolling() {}; - virtual void coeffs_to_local() {}; - virtual void init() {}; virtual void calculate_forces() = 0; }; diff --git a/src/GRANULAR/gran_sub_mod_tangential.h b/src/GRANULAR/gran_sub_mod_tangential.h index b5a195a1d6..27588d49ba 100644 --- a/src/GRANULAR/gran_sub_mod_tangential.h +++ b/src/GRANULAR/gran_sub_mod_tangential.h @@ -63,8 +63,6 @@ class GranSubModTangential : public GranSubMod { public: GranSubModTangential(class GranularModel *, class LAMMPS *); virtual ~GranSubModTangential() {}; - virtual void coeffs_to_local() {}; - virtual void init() {}; virtual void calculate_forces() = 0; double k, damp, mu; // Used by Marshall twisting model }; @@ -105,8 +103,6 @@ class GranSubModTangentialLinearHistoryClassic : public GranSubModTangentialLine public: GranSubModTangentialLinearHistoryClassic(class GranularModel *, class LAMMPS *); void calculate_forces(); - protected: - double xt; }; /* ---------------------------------------------------------------------- */ diff --git a/src/GRANULAR/gran_sub_mod_twisting.h b/src/GRANULAR/gran_sub_mod_twisting.h index c365cba331..46f56d37ce 100644 --- a/src/GRANULAR/gran_sub_mod_twisting.h +++ b/src/GRANULAR/gran_sub_mod_twisting.h @@ -39,8 +39,6 @@ class GranSubModTwisting : public GranSubMod { public: GranSubModTwisting(class GranularModel *, class LAMMPS *); virtual ~GranSubModTwisting() {}; - virtual void coeffs_to_local() {}; - virtual void init() {}; virtual void calculate_forces() = 0; }; @@ -57,8 +55,8 @@ class GranSubModTwistingNone : public GranSubModTwisting { class GranSubModTwistingMarshall : public GranSubModTwisting { public: GranSubModTwistingMarshall(class GranularModel *, class LAMMPS *); + void init() override; void calculate_forces(); - void init(); protected: double k_tang, mu_tang; }; diff --git a/src/GRANULAR/granular_model.cpp b/src/GRANULAR/granular_model.cpp index 5d89e22d03..afb3328d70 100644 --- a/src/GRANULAR/granular_model.cpp +++ b/src/GRANULAR/granular_model.cpp @@ -148,12 +148,13 @@ void GranularModel::construct_submodel(std::string model_name, SubmodelType mode sub_models[model_type]->name.assign(model_name); sub_models[model_type]->allocate_coeffs(); - if (model_type == NORMAL) normal_model = dynamic_cast (sub_models[NORMAL]); - if (model_type == DAMPING) damping_model = dynamic_cast (sub_models[DAMPING]); - if (model_type == TANGENTIAL) tangential_model = dynamic_cast (sub_models[TANGENTIAL]); - if (model_type == ROLLING) rolling_model = dynamic_cast (sub_models[ROLLING]); - if (model_type == TWISTING) twisting_model = dynamic_cast (sub_models[TWISTING]); - if (model_type == HEAT) heat_model = dynamic_cast (sub_models[HEAT]); + // Assign specific submodel pointer + if (model_type == NORMAL) normal_model = dynamic_cast (sub_models[model_type]); + if (model_type == DAMPING) damping_model = dynamic_cast (sub_models[model_type]); + if (model_type == TANGENTIAL) tangential_model = dynamic_cast (sub_models[model_type]); + if (model_type == ROLLING) rolling_model = dynamic_cast (sub_models[model_type]); + if (model_type == TWISTING) twisting_model = dynamic_cast (sub_models[model_type]); + if (model_type == HEAT) heat_model = dynamic_cast (sub_models[model_type]); } /* ---------------------------------------------------------------------- */ @@ -248,7 +249,8 @@ void GranularModel::init() beyond_contact = 1; size_history += sub_models[i]->size_history; if (!sub_models[i]->allow_cohesion && normal_model->cohesive_flag) - error->all(FLERR,"Cannot use {} model with a cohesive normal model, {}", sub_models[i]->name, normal_model->name); + error->all(FLERR,"Cannot use {} model with a cohesive normal model, {}", + sub_models[i]->name, normal_model->name); if (sub_models[i]->area_flag) area_flag = 1; } @@ -260,14 +262,14 @@ void GranularModel::init() int j; for (int i = 0; i < size_history; i++) { - // Find which model owns this history value + // Find which submodel owns this history value size_cumulative = 0; for (j = 0; j < NSUBMODELS; j++) { if ((size_cumulative + sub_models[j]->size_history) > i) break; size_cumulative += sub_models[j]->size_history; } - // Check if model has nondefault transfers, if so copy its array + // Check if submodel has nondefault transfers, if so copy its array transfer_history_factor[i] = -1; if (j != NSUBMODELS) { if (sub_models[j]->nondefault_history_transfer) { @@ -284,8 +286,7 @@ void GranularModel::init() int GranularModel::mix_coeffs(GranularModel *g1, GranularModel *g2) { - int i; - for (i = 0; i < NSUBMODELS; i++) { + for (int i = 0; i < NSUBMODELS; i++) { if (g1->sub_models[i]->name != g2->sub_models[i]->name) return i; construct_submodel(g1->sub_models[i]->name, (SubmodelType) i); @@ -311,6 +312,8 @@ void GranularModel::write_restart(FILE *fp) fwrite(&num_coeffs, sizeof(int), 1, fp); fwrite(sub_models[i]->coeffs, sizeof(double), num_coeffs, fp); } + + fwrite(&limit_damping, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- */ @@ -328,22 +331,23 @@ void GranularModel::read_restart(FILE *fp) if (comm->me == 0) utils::sfread(FLERR, const_cast(model_name.data()), sizeof(char),num_char, fp, nullptr, error); MPI_Bcast(const_cast(model_name.data()), num_char, MPI_CHAR, 0, world); - construct_submodel(model_name, (SubmodelType) i); - if (comm->me == 0) { + if (comm->me == 0) utils::sfread(FLERR, &num_coeff, sizeof(int), 1, fp, nullptr, error); - if (num_coeff != sub_models[i]->num_coeffs) - error->one(FLERR, "Invalid granular model written to restart file"); - } MPI_Bcast(&num_coeff, 1, MPI_INT, 0, world); + if (num_coeff != sub_models[i]->num_coeffs) + error->all(FLERR, "Invalid granular model written to restart file"); - if (comm->me == 0) { + if (comm->me == 0) utils::sfread(FLERR, sub_models[i]->coeffs, sizeof(double), num_coeff, fp, nullptr, error); - } MPI_Bcast(sub_models[i]->coeffs, num_coeff, MPI_DOUBLE, 0, world); sub_models[i]->coeffs_to_local(); } + + if (comm->me == 0) + utils::sfread(FLERR, &limit_damping, sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&limit_damping, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- */ @@ -393,7 +397,7 @@ void GranularModel::calculate_forces() sub3(vi, vj, vr); // normal component - vnnr = dot3(vr, nx); //v_R . n + vnnr = dot3(vr, nx); scale3(vnnr, nx, vn); // tangential component @@ -433,7 +437,7 @@ void GranularModel::calculate_forces() if (contact_type == PAIR) { copy3(torquesi, torquesj); - // Classic pair styles wouldn't scale, but classic option is only used by walls + // Classic pair styles don't scale, but classic option is currently only used by walls dist_to_contact = radi - 0.5 * delta; scale3(dist_to_contact, torquesi); dist_to_contact = radj - 0.5 * delta; diff --git a/src/GRANULAR/granular_model.h b/src/GRANULAR/granular_model.h index de4f82c01e..a047bc27e5 100644 --- a/src/GRANULAR/granular_model.h +++ b/src/GRANULAR/granular_model.h @@ -68,7 +68,7 @@ class GranularModel : protected Pointers { GranSubModRolling *rolling_model; GranSubModTwisting *twisting_model; GranSubModHeat *heat_model; - GranSubMod *sub_models[NSUBMODELS]; // Need to resize if we add more model flavors + GranSubMod *sub_models[NSUBMODELS]; // Extra options int beyond_contact, limit_damping, history_update; @@ -94,9 +94,9 @@ class GranularModel : protected Pointers { bool touch; protected: - int rolling_defined, twisting_defined, heat_defined; // Used to quickly skip undefined submodels - int classic_model; - int area_flag; + int rolling_defined, twisting_defined, heat_defined; // Flag optional submodels + int classic_model; // Flag original pair/gran calculations + int area_flag; // Flag whether area is needed int nclass; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index b6c9bba705..47bfa89da0 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -57,8 +57,6 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) svector = new double[single_extra]; neighprev = 0; - nmodels = 0; - maxmodels = 0; nmax = 0; mass_rigid = nullptr; @@ -302,6 +300,7 @@ void PairGranular::allocate() maxmodels = n * n + 1; // should never need any more space models_list = (GranularModel **) memory->smalloc(maxmodels * sizeof(GranularModel *), "pair:models_list"); for (int i = 0; i < maxmodels; i++) models_list[i] = nullptr; + nmodels = 0; onerad_dynamic = new double[n+1]; onerad_frozen = new double[n+1]; @@ -339,7 +338,6 @@ void PairGranular::coeff(int narg, char **arg) utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); - // Construct new model models_list[nmodels] = new GranularModel(lmp); class GranularModel* model = models_list[nmodels]; @@ -374,9 +372,10 @@ void PairGranular::coeff(int narg, char **arg) } else error->all(FLERR, "Illegal pair_coeff command {}", arg[iarg]); } - // Define default damping model if unspecified, has no coeffs + // Define default damping submodel if unspecified, has no coeffs if (!model->damping_model) model->construct_submodel("viscoelastic", DAMPING); + model->init(); int count = 0; for (int i = ilo; i <= ihi; i++) { @@ -388,7 +387,7 @@ void PairGranular::coeff(int narg, char **arg) } } - // If there are > ntype^2 models, delete unused ones + // If there are > ntype^2 models, delete unused models if (nmodels == maxmodels) prune_models(); if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); @@ -419,7 +418,6 @@ void PairGranular::init_style() int size_max[NSUBMODELS] = {0}; for (int n = 0; n < nmodels; n++) { model = models_list[n]; - model->init(); if (model->beyond_contact) { beyond_contact = 1; @@ -559,7 +557,6 @@ double PairGranular::init_one(int i, int j) model1->sub_models[error_code]->name, model2->sub_models[error_code]->name); - // Initialize new model, init_one() occurs after init_style model->init(); for (int k = 0; k < NSUBMODELS; k++) @@ -657,8 +654,10 @@ void PairGranular::read_restart(FILE *fp) MPI_Bcast(&nmodels,1,MPI_INT,0,world); for (i = 0; i < nmodels; i++) { + delete models_list[i]; models_list[i] = new GranularModel(lmp); models_list[i]->read_restart(fp); + models_list[i]->init(); } for (i = 1; i <= atom->ntypes; i++) { @@ -690,28 +689,22 @@ double PairGranular::single(int i, int j, int itype, int jtype, double rsq, double /* factor_coul */, double factor_lj, double &fforce) { - if (factor_lj == 0) { fforce = 0.0; for (int m = 0; m < single_extra; m++) svector[m] = 0.0; return 0.0; } - int jnum; - int *jlist; - double *history,*allhistory; - int nall = atom->nlocal + atom->nghost; if ((i >= nall) || (j >= nall)) error->all(FLERR,"Not enough atoms for pair granular single function"); - int touchflag; - double **x = atom->x; - double *radius = atom->radius; - class GranularModel* model = models_list[types_indices[itype][jtype]]; // Reset model and copy initial geometric data + double **x = atom->x; + double *radius = atom->radius; + model->xi = x[i]; model->xj = x[j]; model->radi = radius[i]; @@ -719,8 +712,9 @@ double PairGranular::single(int i, int j, int itype, int jtype, model->history_update = 0; // Don't update history // If history is needed - jnum = list->numneigh[i]; - jlist = list->firstneigh[i]; + double *history,*allhistory; + int jnum = list->numneigh[i]; + int *jlist = list->firstneigh[i]; if (use_history) { if ((fix_history == nullptr) || (fix_history->firstvalue == nullptr)) error->one(FLERR,"Pair granular single computation needs history"); @@ -730,11 +724,11 @@ double PairGranular::single(int i, int j, int itype, int jtype, if (neighprev >= jnum) neighprev = 0; if (jlist[neighprev] == j) break; } - history = &allhistory[size_history*neighprev]; + history = &allhistory[size_history * neighprev]; model->touch = fix_history->firstflag[i][neighprev]; } - touchflag = model->check_contact(); + int touchflag = model->check_contact(); if (!touchflag) { fforce = 0.0; @@ -745,17 +739,16 @@ double PairGranular::single(int i, int j, int itype, int jtype, // meff = effective mass of pair of particles // if I or J part of rigid body, use body mass // if I or J is frozen, meff is other particle - double mi, mj, meff; double *rmass = atom->rmass; int *mask = atom->mask; - mi = rmass[i]; - mj = rmass[j]; + double mi = rmass[i]; + double mj = rmass[j]; if (fix_rigid) { if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; } - meff = mi * mj / (mi + mj); + double meff = mi * mj / (mi + mj); if (mask[i] & freeze_group_bit) meff = mj; if (mask[j] & freeze_group_bit) meff = mi; @@ -770,11 +763,10 @@ double PairGranular::single(int i, int j, int itype, int jtype, model->omegaj = omega[j]; model->history = history; - double *forces, *torquesi, *torquesj; model->calculate_forces(); - forces = model->forces; - torquesi = model->torquesi; - torquesj = model->torquesj; + double *forces = model->forces; + double *torquesi = model->torquesi; + double *torquesj = model->torquesj; // apply forces & torques fforce = MathExtra::len3(forces); @@ -843,7 +835,7 @@ void PairGranular::transfer_history(double* source, double* target, int itype, i class GranularModel* model = models_list[types_indices[itype][jtype]]; if (model->nondefault_history_transfer) { for (int i = 0; i < size_history; i++) { - target[i] = model->transfer_history_factor[i] * source [i]; + target[i] = model->transfer_history_factor[i] * source[i]; } } else { for (int i = 0; i < size_history; i++) {