diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index f8d351051c..16b1b228ab 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -87,7 +87,7 @@ class FixWallGran : public Fix { int store; - + }; } // namespace LAMMPS_NS diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index 7eba3bdfb3..ba8cd467ec 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/src/GRANULAR/fix_wall_gran_region.cpp @@ -231,7 +231,7 @@ void FixWallGranRegion::post_force(int /*vflag*/) model->radi = radius[i]; model->radj = region->contact[ic].radius; model->r = region->contact[ic].r; - model->i = i; // Added for MDR + model->i = i; // Added for MDR model->j = ic; // Added for MDR if (model->beyond_contact) model->touch = history_many[i][c2r[ic]][0]; diff --git a/src/GRANULAR/fix_wall_gran_region.h b/src/GRANULAR/fix_wall_gran_region.h index e9e24415bc..6f9866d8ae 100644 --- a/src/GRANULAR/fix_wall_gran_region.h +++ b/src/GRANULAR/fix_wall_gran_region.h @@ -53,20 +53,20 @@ class FixWallGranRegion : public FixWallGran { int tmax; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL private: - + int nregion; // shear history for multiple contacts per particle - - - + + + // c2r[i] = index of Ith contact in // region-contact[] list of contacts int motion_resetflag; // used by restart to indicate that region // vel info is to be reset - + }; } // namespace LAMMPS_NS diff --git a/src/GRANULAR/gran_sub_mod_normal.cpp b/src/GRANULAR/gran_sub_mod_normal.cpp index bf3d045c2f..3521f0be72 100644 --- a/src/GRANULAR/gran_sub_mod_normal.cpp +++ b/src/GRANULAR/gran_sub_mod_normal.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include using namespace LAMMPS_NS; @@ -401,12 +401,12 @@ 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; + size_history = 26; nondefault_history_transfer = 1; transfer_history_factor = new double[size_history]; //transfer_history_factor[0] = +1; - for (int i = 0; i < size_history; i++) { + for (int i = 0; i < size_history; i++) { transfer_history_factor[i] = +1; } } @@ -435,18 +435,18 @@ void GranSubModNormalMDR::coeffs_to_local() double GranSubModNormalMDR::calculate_forces() { - // To understand the structure of the overall code it is important to consider + // To understand the structure of the overall code it is important to consider // the following: // - // The MDR contact model was developed by imagining individual particles being - // squished between a number of rigid flats (references below). To allow - // for many interacting particles, we extend the idea of isolated particles surrounded - // by rigid flats. In particular, we imagine placing rigid flats at the overlap + // The MDR contact model was developed by imagining individual particles being + // squished between a number of rigid flats (references below). To allow + // for many interacting particles, we extend the idea of isolated particles surrounded + // by rigid flats. In particular, we imagine placing rigid flats at the overlap // midpoints between particles. The force is calculated seperately on both sides - // of the contact assuming interaction with a rigid flat. The two forces are then - // averaged on either side of the contact to determine the final force. If the + // of the contact assuming interaction with a rigid flat. The two forces are then + // averaged on either side of the contact to determine the final force. If the // contact is between a particle and wall then only one force evaluation is required. - // + // // Zunker and Kamrin, 2024, Part I: https://doi.org/10.1016/j.jmps.2023.105492 // Zunker and Kamrin, 2024, Part II: https://doi.org/10.1016/j.jmps.2023.105493 @@ -456,8 +456,8 @@ double GranSubModNormalMDR::calculate_forces() const int j_true = gm->j; // true j particle index const double radi_true = gm->radi; // true i particle initial radius const double radj_true = gm->radj; // true j particle initial radius - - F = 0.0; // average force + + F = 0.0; // average force double F0 = 0.0; // force on contact side 0 double F1 = 0.0; // force on contact side 1 double R0 = 0.0; @@ -465,21 +465,21 @@ double GranSubModNormalMDR::calculate_forces() int i0 = 0; int i1 = 0; double delta = gm->delta; // apparent overlap - - //if (gm->contact_type == PAIR) delta = gm->delta/2.0; // half displacement to imagine interaction with rigid flat - //std::cout << "Normal force is called for: " << i_true << ", " << j_true << std::endl; - //std::cout << "Contact model has been entered " << gm->contact_type << ", " << PAIR << ", " << WALL << ", " << WALLREGION << ", " << gm->itype << ", " << gm->jtype << ", " << gm->delta << std::endl; - // initialize indexing in history array of different constact history variables - const int delta_offset_0 = 0; // apparent overlap - const int delta_offset_1 = 1; - const int deltao_offset_0 = 2; // displacement - const int deltao_offset_1 = 3; + //if (gm->contact_type == PAIR) delta = gm->delta/2.0; // half displacement to imagine interaction with rigid flat + //std::cout << "Normal force is called for: " << i_true << ", " << j_true << std::endl; + //std::cout << "Contact model has been entered " << gm->contact_type << ", " << PAIR << ", " << WALL << ", " << WALLREGION << ", " << gm->itype << ", " << gm->jtype << ", " << gm->delta << std::endl; + + // initialize indexing in history array of different constact history variables + const int delta_offset_0 = 0; // apparent overlap + const int delta_offset_1 = 1; + const int deltao_offset_0 = 2; // displacement + const int deltao_offset_1 = 3; const int delta_MDR_offset_0 = 4; // MDR apparent overlap const int delta_MDR_offset_1 = 5; const int delta_BULK_offset_0 = 6; // bulk displacement const int delta_BULK_offset_1 = 7; - const int deltamax_MDR_offset_0 = 8; // maximum MDR apparent overlap + const int deltamax_MDR_offset_0 = 8; // maximum MDR apparent overlap const int deltamax_MDR_offset_1 = 9; const int Yflag_offset_0 = 10; // yield flag const int Yflag_offset_1 = 11; @@ -493,32 +493,32 @@ double GranSubModNormalMDR::calculate_forces() const int Ac_offset_1 = 19; const int eps_bar_offset_0 = 20; // volume-averaged infinitesimal strain tensor const int eps_bar_offset_1 = 21; - const int penalty_offset_ = 22; // contact penalty + const int penalty_offset_ = 22; // contact penalty const int deltamax_offset_ = 23; const int deltap_offset_0 = 24; const int deltap_offset_1 = 25; - // initialize particle history variables + // 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_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_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 = 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 + 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]; @@ -527,7 +527,7 @@ double GranSubModNormalMDR::calculate_forces() 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 * 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]; @@ -537,9 +537,9 @@ double GranSubModNormalMDR::calculate_forces() 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 + + double * history = & gm->history[history_index]; // load in all history variables // Rigid flat placement scheme double * deltamax_offset = & history[deltamax_offset_]; @@ -552,19 +552,19 @@ double GranSubModNormalMDR::calculate_forces() if (gm->delta >= *deltamax_offset) *deltamax_offset = gm->delta; deltamax = *deltamax_offset; - for (int contactSide = 0; contactSide < 2; contactSide++) { + for (int contactSide = 0; contactSide < 2; contactSide++) { - double * delta_offset; + double * delta_offset; double * deltao_offset; - double * delta_MDR_offset; - double * delta_BULK_offset; - double * deltamax_MDR_offset; - double * Yflag_offset; - double * deltaY_offset; + double * delta_MDR_offset; + double * delta_BULK_offset; + double * deltamax_MDR_offset; + double * Yflag_offset; + double * deltaY_offset; double * cA_offset; - double * aAdh_offset; - double * Ac_offset; - double * eps_bar_offset; + double * aAdh_offset; + double * Ac_offset; + double * eps_bar_offset; double * penalty_offset; double * deltap_offset; @@ -631,7 +631,7 @@ double GranSubModNormalMDR::calculate_forces() gm->radi = radj_true; gm->radj = radi_true; } - + double delta_geo; double delta_geo_alt; double delta_geoOpt1 = deltamax*(deltamax - 2.0*R1)/(2.0*(deltamax - R0 - R1)); @@ -649,7 +649,7 @@ double GranSubModNormalMDR::calculate_forces() delta = delta_geo + (deltap1 - delta_geo)/(deltap - deltamax)*(gm->delta-deltamax); //std::cout << "CS 1: " << gm->radi << ", " << gm->radj << ", gm->delta " << gm->delta << ", delta " << delta << ", deltamax " << deltamax << ", delta_geo " << delta_geo << ", delta_geo_alt " << delta_geo_alt << ", delta_geo/Ri " << delta_geo/gm->radi << ", delta_geo_alt/Rj " << delta_geo_alt/gm->radj << ", delta_geoOpt1 " << delta_geoOpt1 << ", delta_geoOpt2 " << delta_geoOpt2 << ", deltamax-delta " << (deltamax-gm->delta) << std::endl; - + delta_offset = & history[delta_offset_1]; deltao_offset = & history[deltao_offset_1]; delta_MDR_offset = & history[delta_MDR_offset_1]; @@ -678,8 +678,8 @@ double GranSubModNormalMDR::calculate_forces() const double Ro = Rinitial[i]; // initial radius const double R = gm->radi; // apparent radius - - // kinematics + + // kinematics const double ddelta = delta - *delta_offset; *delta_offset = delta; @@ -689,7 +689,7 @@ double GranSubModNormalMDR::calculate_forces() double ddelta_MDR; double ddelta_BULK; - if ( psi[i] < psi_b ) { // if true, bulk response has triggered, split displacement increment between the MDR and BULK components + if ( psi[i] < psi_b ) { // if true, bulk response has triggered, split displacement increment between the MDR and BULK components ddelta_MDR = std::min(ddelta-ddelta_bar[i], delta-*delta_MDR_offset); ddelta_BULK = ddelta_bar[i]; } else { // if false, no bulk response, full displacement increment goes to the MDR component @@ -712,7 +712,7 @@ double GranSubModNormalMDR::calculate_forces() *deltaY_offset = delta_MDR; *cA_offset = M_PI*(pow(*deltaY_offset,2.0) - *deltaY_offset*R); } - } + } //if (i_true == 167 && j_true == 204) { //std::cout << "i " << i << " | j " << j << " | delta_BULK: " << delta_BULK << " | delta_MDR " << delta_MDR << " | ddelta_BULK " << ddelta_BULK << " | ddelta_MDR " << ddelta_MDR << std::endl; @@ -725,25 +725,25 @@ double GranSubModNormalMDR::calculate_forces() double A; // height of elliptical indenter double B; // width of elliptical indenter double deltae1D; // transformed elastic displacement - double deltaR; // displacement correction + double deltaR; // displacement correction double amax; // maximum experienced contact radius const double cA = *cA_offset; // contact area intercept if ( *Yflag_offset == 0.0 ) { // elastic contact - A = 4.0*R; - B = 2.0*R; - deltae1D = delta_MDR; - (deltae1D > 0) ? amax = sqrt(deltae1D*R) : amax = 0.0; + A = 4.0*R; + B = 2.0*R; + deltae1D = delta_MDR; + (deltae1D > 0) ? amax = sqrt(deltae1D*R) : amax = 0.0; } else { // plastic contact - amax = sqrt((2.0*deltamax_MDR*R - pow(deltamax_MDR,2.0)) + cA/M_PI); - A = 4.0*pY/Eeff*amax; - B = 2.0*amax; - const double deltae1Dmax = A/2.0; // maximum transformed elastic displacement + amax = sqrt((2.0*deltamax_MDR*R - pow(deltamax_MDR,2.0)) + cA/M_PI); + A = 4.0*pY/Eeff*amax; + B = 2.0*amax; + const double deltae1Dmax = A/2.0; // maximum transformed elastic displacement const double Fmax = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltae1Dmax/A) - (1.0 - 2.0*deltae1Dmax/A)*sqrt(4.0*deltae1Dmax/A - 4.0*pow(deltae1Dmax,2.0)/pow(A,2.0))); // force caused by full submersion of elliptical indenter to depth of A/2 const double zR = R - (deltamax_MDR - deltae1Dmax); // depth of particle center - deltaR = (Fmax*(2*pow(amax,2.0)*(-1 + nu) - (-1 + 2*nu)*zR*(-zR + sqrt(pow(amax,2.0) + pow(zR,2.0)))))/((M_PI*pow(amax,2.0))*2*G*sqrt(pow(amax,2.0) + pow(zR,2.0))); - deltae1D = (delta_MDR - deltamax_MDR + deltae1Dmax + deltaR)/(1 + deltaR/deltae1Dmax); // transformed elastic displacement - + deltaR = (Fmax*(2*pow(amax,2.0)*(-1 + nu) - (-1 + 2*nu)*zR*(-zR + sqrt(pow(amax,2.0) + pow(zR,2.0)))))/((M_PI*pow(amax,2.0))*2*G*sqrt(pow(amax,2.0) + pow(zR,2.0))); + deltae1D = (delta_MDR - deltamax_MDR + deltae1Dmax + deltaR)/(1 + deltaR/deltae1Dmax); // transformed elastic displacement + // added for rigid flat placement *deltap_offset = deltamax_MDR - (deltae1Dmax + deltaR); //std::cout << *deltap_offset << std::endl; @@ -755,7 +755,7 @@ double GranSubModNormalMDR::calculate_forces() double a_na; double a_fac = 0.99; (deltae1D >= 0.0) ? a_na = B*sqrt(A - deltae1D)*sqrt(deltae1D)/A : a_na = 0.0; - double aAdh = *aAdh_offset; + double aAdh = *aAdh_offset; if (aAdh > a_fac*amax) aAdh = a_fac*amax; //if (i_true == 4 && j_true == 52){ @@ -764,9 +764,9 @@ double GranSubModNormalMDR::calculate_forces() if ( gamma > 0.0 ) { // adhesive contact double g_aAdh; - + if (delta_MDR == deltamax_MDR || a_na >= aAdh ) { // case 1: no tensile springs, purely compressive contact - (deltae1D <= 0.0) ? F_MDR = 0.0 : F_MDR = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltae1D/A) - (1.0 - 2.0*deltae1D/A)*sqrt(4.0*deltae1D/A - 4.0*pow(deltae1D,2.0)/pow(A,2.0))); + (deltae1D <= 0.0) ? F_MDR = 0.0 : F_MDR = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltae1D/A) - (1.0 - 2.0*deltae1D/A)*sqrt(4.0*deltae1D/A - 4.0*pow(deltae1D,2.0)/pow(A,2.0))); if ( std::isnan(F_MDR) ) { std::cout << "F_MDR is NaN, case 1: no tensile springs" << std::endl; //std::cout << "Normal model: " << gm->delta << ", " << ddelta << ", " << gm->radi << ", " << gm->radj << " | delta: " << delta0 << ", " << delta1 << " | delta2_offset: " << *delta2_offset0 << ", " << *delta2_offset1 << "| dde: " << dde0 << ", " << dde1 << "| Fold: " << F0old << ", " << F1old << " | a: " << a0 << ", " << a1 << " | k_BULK: " << k_BULK0 << ", " << k_BULK1 << " | h_BULK: " << h_BULK0 << ", " << h_BULK1 << std::endl; @@ -775,7 +775,7 @@ double GranSubModNormalMDR::calculate_forces() } *aAdh_offset = a_fac*a_na; } else { - const double lmax = sqrt(2.0*M_PI*aAdh*gamma/Eeff); + const double lmax = sqrt(2.0*M_PI*aAdh*gamma/Eeff); g_aAdh = A/2 - A/B*sqrt(pow(B,2.0)/4 - pow(aAdh,2.0)); const double acrit = (-((pow(B,2)*gamma*M_PI)/(pow(A,2)*Eeff)) + (pow(2,0.3333333333333333)*pow(B,4)*pow(gamma,2)*pow(M_PI,1.6666666666666667))/ (pow(A,2)*pow(Eeff,2)*pow((27*pow(A,4)*pow(B,4)*gamma)/Eeff - (2*pow(B,6)*pow(gamma,3)*pow(M_PI,2))/pow(Eeff,3) + (3*sqrt(3)*sqrt(27*pow(A,8)*pow(B,8)*pow(Eeff,2)*pow(gamma,2) - @@ -785,19 +785,19 @@ double GranSubModNormalMDR::calculate_forces() if ( (deltae1D + lmax - g_aAdh) >= 0.0) { // case 2: tensile springs, but not exceeding critical length --> deltae + lmax - g(aAdhes) >= 0 //if (contactSide == 0) { - //std::cout << "Case 2 tensile springs not exceeding critical length, R " << R << "deltae1D " << deltae1D << " , lmax " << lmax << " , g_adh" << g_aAdh << " sum " << (deltae1D + lmax - g_aAdh) << std::endl; + //std::cout << "Case 2 tensile springs not exceeding critical length, R " << R << "deltae1D " << deltae1D << " , lmax " << lmax << " , g_adh" << g_aAdh << " sum " << (deltae1D + lmax - g_aAdh) << std::endl; //} - const double deltaeAdh = g_aAdh; + const double deltaeAdh = g_aAdh; const double F_na = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltaeAdh/A) - (1.0 - 2.0*deltaeAdh/A)*sqrt(4.0*deltaeAdh/A - 4.0*pow(deltaeAdh,2.0)/pow(A,2.0))); const double F_Adhes = 2.0*Eeff*(deltae1D - deltaeAdh)*aAdh; - F_MDR = F_na + F_Adhes; + F_MDR = F_na + F_Adhes; if ( std::isnan(F_MDR) ) std::cout << "F_MDR is NaN, case 2: tensile springs, but not exceeding critical length" << std::endl; } else { // case 3: tensile springs exceed critical length --> deltae + lmax - g(aAdhes) = 0 //if (contactSide == 0) { - //std::cout << "Case 3 tensile springs exceed critical length, R " << R << " deltae1D " << deltae1D << " , lmax " << lmax << " , g_adh " << g_aAdh << " sum " << (deltae1D + lmax - g_aAdh) << std::endl; + //std::cout << "Case 3 tensile springs exceed critical length, R " << R << " deltae1D " << deltae1D << " , lmax " << lmax << " , g_adh " << g_aAdh << " sum " << (deltae1D + lmax - g_aAdh) << std::endl; //} - + if ( aAdh < acrit ) { aAdh = 0.0; F_MDR = 0.0; @@ -807,31 +807,31 @@ double GranSubModNormalMDR::calculate_forces() const double error = 1e-10; const double error2 = 1e-16; double aAdh_tmp = aAdh; - double fa; + double fa; double fa2; double dfda; for (int lv1 = 0; lv1 < maxIterations; ++lv1) { fa = deltae1D + sqrt(2.0*M_PI*aAdh_tmp*gamma/Eeff) - ( A/2 - A/B*sqrt(pow(B,2.0)/4 - pow(aAdh_tmp,2.0)) ); if (abs(fa) < error) { break; - } + } dfda = -((aAdh_tmp*A)/(B*sqrt(-pow(aAdh_tmp,2.0) + pow(B,2.0)/4.0))) + (gamma*sqrt(M_PI/2.0))/(Eeff*sqrt((aAdh_tmp*gamma)/Eeff)); aAdh_tmp = aAdh_tmp - fa/dfda; fa2 = deltae1D + sqrt(2.0*M_PI*aAdh_tmp*gamma/Eeff) - ( A/2 - A/B*sqrt(pow(B,2.0)/4 - pow(aAdh_tmp,2.0)) ); if (abs(fa-fa2) < error2) { break; - } + } if (lv1 == maxIterations-1){ aAdh_tmp = 0.0; } } aAdh = aAdh_tmp; - - g_aAdh = A/2.0 - A/B*sqrt(pow(B,2.0)/4.0 - pow(aAdh,2.0)); - const double deltaeAdh = g_aAdh; + + g_aAdh = A/2.0 - A/B*sqrt(pow(B,2.0)/4.0 - pow(aAdh,2.0)); + const double deltaeAdh = g_aAdh; const double F_na = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltaeAdh/A) - (1.0 - 2.0*deltaeAdh/A)*sqrt(4.0*deltaeAdh/A - 4.0*pow(deltaeAdh,2.0)/pow(A,2.0))); const double F_Adhes = 2.0*Eeff*(deltae1D - deltaeAdh)*aAdh; - F_MDR = F_na + F_Adhes; + F_MDR = F_na + F_Adhes; if ( std::isnan(F_MDR) ) std::cout << "F_MDR is NaN, case 3: tensile springs exceed critical length" << std::endl; } *aAdh_offset = aAdh; @@ -839,13 +839,13 @@ double GranSubModNormalMDR::calculate_forces() } } else { // non-adhesive contact *aAdh_offset = a_na; - (deltae1D <= 0.0) ? F_MDR = 0.0 : F_MDR = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltae1D/A) - (1.0 - 2.0*deltae1D/A)*sqrt(4.0*deltae1D/A - 4.0*pow(deltae1D,2.0)/pow(A,2.0))); + (deltae1D <= 0.0) ? F_MDR = 0.0 : F_MDR = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltae1D/A) - (1.0 - 2.0*deltae1D/A)*sqrt(4.0*deltae1D/A - 4.0*pow(deltae1D,2.0)/pow(A,2.0))); if ( std::isnan(F_MDR) ) { std::cout << "F_MDR is NaN, non-adhesive case" << std::endl; //std::cout << "Normal model: " << gm->delta << ", " << ddelta << ", " << gm->radi << ", " << gm->radj << " | delta: " << delta0 << ", " << delta1 << " | delta2_offset: " << *delta2_offset0 << ", " << *delta2_offset1 << "| dde: " << dde0 << ", " << dde1 << "| Fold: " << F0old << ", " << F1old << " | a: " << a0 << ", " << a1 << " | k_BULK: " << k_BULK0 << ", " << k_BULK1 << " | h_BULK: " << h_BULK0 << ", " << h_BULK1 << std::endl; std::cout << "i_true: " << i_true << ", j_true: " << j_true << ", i_tag: " << atom->tag[i_true] << ", j_tag: " << atom->tag[j_true] << ", deltae1D: " << deltae1D << ", A: " << A << ", B: " << B << ", amax: " << amax << ", deltamax_MDR: " << deltamax_MDR << ", R: " << R << std::endl; std::exit(1); - } + } } //std::cout << gm->i << ", " << gm->j << ", aAdh_offset: " << *aAdh_offset << ", aAdh: " << aAdh << ", a_na: " << a_na << std::endl; @@ -860,8 +860,8 @@ double GranSubModNormalMDR::calculate_forces() //std::cout << gm->i << ", " << gm->j << ", " << gm->contact_type << ", " << *gm->xi[1] << ", " << *gm->xj[1] << std::endl; - // area related calculations - double Ac; + // area related calculations + double Ac; (*Yflag_offset == 0.0) ? Ac = M_PI*delta*R : Ac = M_PI*((2.0*delta*R - pow(delta,2.0)) + cA/M_PI); if (Ac < 0.0 ) Ac = 0.0; Atot_sum[i] += wij*(Ac - 2.0*M_PI*R*(deltamax_MDR + delta_BULK)); @@ -932,7 +932,7 @@ double GranSubModNormalMDR::calculate_forces() //int rank = 0; //MPI_Comm_rank(MPI_COMM_WORLD, &rank); - //std::cout << "Step: " << lmp->update->ntimestep << ", CS: " << contactSide << ", itag: " << atom->tag[i] << ", jtag: " << atom->tag[j] << ", rank: " << rank << ", gm->delta: " << gm->delta << ", delta: " << delta << ", ddelta: " << ddelta << ", delta_bar: " << ddelta_bar[i] << ", R: " << R << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << std::endl; + //std::cout << "Step: " << lmp->update->ntimestep << ", CS: " << contactSide << ", itag: " << atom->tag[i] << ", jtag: " << atom->tag[j] << ", rank: " << rank << ", gm->delta: " << gm->delta << ", delta: " << delta << ", ddelta: " << ddelta << ", delta_bar: " << ddelta_bar[i] << ", R: " << R << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << std::endl; //std::cout << gm->i << ", " << gm->j << ", " << Vgeo[i] << ", " << Acon0[i] << ", " << Acon1[i] << ", " << Ac << ", " << kappa << " || " << psi[i] << ", " << ddelta_bar[i] << ", " << ddelta << ", " << ddelta_MDR << ", " << ddelta_BULK << ", " << delta << ", " << delta_MDR << ", " << delta_BULK << ", " << F_MDR << ", " << F_BULK << ", " << R << " || " << deltae1D << ", " << A << ", " << B << std::endl; @@ -950,7 +950,7 @@ double GranSubModNormalMDR::calculate_forces() // radius update scheme quantity calculation Vcaps[i] += (M_PI/3.0)*pow(delta,2.0)*(3.0*R - delta); - + const double Fntmp = wij*(F_MDR + F_BULK); const double fx = Fntmp*gm->nx[0]; const double fy = Fntmp*gm->nx[1]; @@ -960,7 +960,7 @@ double GranSubModNormalMDR::calculate_forces() const double bz = -(Ro - deltao)*gm->nx[2]; const double eps_bar_contact = (1.0/(3*kappa*Velas[i]))*(fx*bx + fy*by + fz*bz); eps_bar[i] += eps_bar_contact; - + //double **x = atom->x; //const double xi = x[i_true][0]; //const double xj = x[j_true][0]; @@ -977,7 +977,7 @@ double GranSubModNormalMDR::calculate_forces() //std::cout << i << ", " << j << ", " << gm->contact_type << " || " << eps_bar_contact << ", " << *eps_bar_offset << ", " << (uintptr_t)(eps_bar_offset) << " || " << wij << ", " << ddeltao << ", " << deltao << " || " << Ro << ", " << R << std::endl; //} - + double desp_bar_contact = eps_bar_contact - *eps_bar_offset; // && desp_bar_contact < 0.0 if(delta_MDR == deltamax_MDR && *Yflag_offset > 0.0 && F_MDR > 0.0){ const double Vo = (4.0/3.0)*M_PI*pow(Ro,3.0); @@ -1020,7 +1020,7 @@ double GranSubModNormalMDR::calculate_forces() //} } - + gm->i = i_true; gm->j = j_true; gm->radi = radi_true; @@ -1044,7 +1044,7 @@ double GranSubModNormalMDR::calculate_forces() if (gm->contact_type != PAIR) { F = wij*F0*wallForceMagnifer; } else { - F = wij*(F0 + F1)/2.0; + F = wij*(F0 + F1)/2.0; } //std::cout << F << ", " << F0 << ", " << F1 << " | " << R0 << ", " << R1 << std::endl; @@ -1073,7 +1073,7 @@ double GranSubModNormalMDR::calculate_forces() //const double xi = x[gm->i][0]; //const double xj = x[gm->j][0]; //const double del = 20.0 - abs(xi-xj); - + //if (i_true == 146 && j_true == 152) { // CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/avicelTableting/rigid_flat_output.csv"); // std::stringstream rowDataStream; @@ -1133,21 +1133,21 @@ void GranSubModNormalMDR::set_fncrit() //std::cout << "F_MDR should be > 0: " << F_MDR << ", " << gm->i << ", " << gm->j << std::endl; //std::cout << gm->i << ", " << gm->j << " || " << delta << ", " << delta_MDR << ", " << deltamax_MDR << ", " << deltae1D << " || " << A << ", " << B << std::endl; -//std::cout << i_true << ", " << j_true << std::endl; +//std::cout << i_true << ", " << j_true << std::endl; //std::cout << history_index << ", " << history[0] << ", " << history[1] << ", " << history[2] << std::endl; // initialize all history variables - //double delta_offset; + //double delta_offset; //double deltao_offset; - //double delta_MDR_offset; - //double delta_BULK_offset; - //double deltamax_MDR_offset; - //double Yflag; - //double deltaY_offset; - //double Ac_offset; - //double aAdh_offset; - //double deltap_offset; - //double cA_offset; + //double delta_MDR_offset; + //double delta_BULK_offset; + //double deltamax_MDR_offset; + //double Yflag; + //double deltaY_offset; + //double Ac_offset; + //double aAdh_offset; + //double deltap_offset; + //double cA_offset; //double eps_bar_offset; - //double wall_contact_flag_offset; \ No newline at end of file + //double wall_contact_flag_offset; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 04d06c296f..230cba446b 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -504,7 +504,7 @@ void PairGranular::init_style() index_vq = atom->find_custom("viscous_heat",tmp1,tmp2); index_vt = atom->find_custom("viscous_temp",tmp1,tmp2); fix_flag = 1; - } + } */ if (model->normal_model->name == "mdr") { @@ -521,22 +521,22 @@ void PairGranular::init_style() 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_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_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 = 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_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_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 contacts on particle std::cout << "MDR history variables have been initialized 2" << ", " << index_Ro << std::endl; @@ -548,7 +548,7 @@ void PairGranular::init_style() modify->add_fix("fix_mdr_mean_surf_disp all mdr/mean/surf/disp"); fix_flag = 1; - } + } } // check for FixFreeze and set freeze_group_bit diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 16a97250b8..10b581968a 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -74,7 +74,7 @@ class PairGranular : public Pair { int index_Ro; int index_Vcaps; int index_Vgeo; - int index_Velas; + int index_Velas; int index_eps_bar; int index_dRnumerator; int index_dRdenominator; diff --git a/src/csv_writer.h b/src/csv_writer.h index 0a01604393..29bb44e419 100644 --- a/src/csv_writer.h +++ b/src/csv_writer.h @@ -22,4 +22,4 @@ public: private: std::string filename_; -}; \ No newline at end of file +}; diff --git a/src/fix_mdr_mean_surf_disp.cpp b/src/fix_mdr_mean_surf_disp.cpp index 4165299b05..cdc170f968 100644 --- a/src/fix_mdr_mean_surf_disp.cpp +++ b/src/fix_mdr_mean_surf_disp.cpp @@ -63,9 +63,9 @@ int FixMDRmeanSurfDisp::setmask() void FixMDRmeanSurfDisp::setup_pre_force(int /*vflag*/) { int tmp1, tmp2; - int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); - int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); - Acon0 = atom->dvector[index_Acon0]; + int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); + int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); + Acon0 = atom->dvector[index_Acon0]; ddelta_bar = atom->dvector[index_ddelta_bar]; pre_force(0); @@ -97,9 +97,9 @@ void FixMDRmeanSurfDisp::pre_force(int) FixNeighHistory * fix_history = dynamic_cast(modify->get_fix_by_id("NEIGH_HISTORY_GRANULAR")); PairGranular * pair = dynamic_cast(force->pair_match("granular",1)); NeighList * list = pair->list; - + const int size_history = pair->get_size_history(); - + //std::cout << " " << std::endl; //std::cout << "New Step" << std::endl; @@ -135,14 +135,14 @@ void FixMDRmeanSurfDisp::pre_force(int) const double ytmp = x[i][1]; const double ztmp = x[i][2]; allhistory = firsthistory[i]; - double radi = radius[i]; + 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]; + 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; @@ -153,7 +153,7 @@ void FixMDRmeanSurfDisp::pre_force(int) const double deltan_ij = radsum_ij - r_ij; if (deltan_ij >= 0.0) { for (int kk = jj; kk < jnum; kk++) { - k = jlist[kk]; + k = jlist[kk]; k &= NEIGHMASK; ktype = type[k]; if (kk != jj) { @@ -175,7 +175,7 @@ void FixMDRmeanSurfDisp::pre_force(int) 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 @@ -186,7 +186,7 @@ void FixMDRmeanSurfDisp::pre_force(int) // 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; + double * pjk = NULL; int * const jklist = firstneigh[j]; const int jknum = numneigh[j]; for (int jk = 0; jk < jknum; jk++) { @@ -219,7 +219,7 @@ void FixMDRmeanSurfDisp::pre_force(int) //std::cout << "Print 198 pjk[0]: " << rank << ", " << pjk[0] << std::endl; break; } - } + } } //std::cout << "Print: " << __LINE__ << std::endl; @@ -234,8 +234,8 @@ void FixMDRmeanSurfDisp::pre_force(int) 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/M_PI - 1.0/2.0)) ); + 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/M_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; @@ -243,7 +243,7 @@ void FixMDRmeanSurfDisp::pre_force(int) 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); + 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/M_PI - 1.0/2.0)) ); } else { // the central particle is i if (j < atom->nlocal || k < atom->nlocal) { @@ -253,7 +253,7 @@ void FixMDRmeanSurfDisp::pre_force(int) 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); + const double alpha = std::acos(enx_ij*enx_ik + eny_ij*eny_ik + enz_ij*enz_ik); //std::cout << "Print: " << __LINE__ << std::endl; pjk[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/M_PI - 1.0/2.0)) ); //std::cout << "Print: " << __LINE__ << std::endl; @@ -264,7 +264,7 @@ void FixMDRmeanSurfDisp::pre_force(int) } } } - } + } } @@ -286,7 +286,7 @@ void FixMDRmeanSurfDisp::pre_force(int) double *radius = atom->radius; int nlocal = atom->nlocal; - + inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -336,8 +336,8 @@ void FixMDRmeanSurfDisp::pre_force(int) //const double wij = std::max(1.0-pij,0.0); const double delta = model->radsum - sqrt(model->rsq); - const int delta_offset_0 = 0; // apparent overlap - const int delta_offset_1 = 1; + const int delta_offset_0 = 0; // apparent overlap + const int delta_offset_1 = 1; const int Ac_offset_0 = 18; // contact area const int Ac_offset_1 = 19; const int deltamax_offset_ = 23; @@ -347,7 +347,7 @@ void FixMDRmeanSurfDisp::pre_force(int) double deltamax = history[deltamax_offset_]; double deltap0 = history[deltap_offset_0]; double deltap1 = history[deltap_offset_1]; - + if (delta > deltamax) deltamax = delta; double delta0old = history[delta_offset_0]; @@ -511,13 +511,13 @@ void FixMDRmeanSurfDisp::pre_force(int) if (Acon0[i] != 0.0) { const double delta = model->radsum - model->r; const double delta_offset0 = fix->history_many[i][fix->c2r[ic]][0]; - const double ddelta = delta - delta_offset0; + const double ddelta = delta - delta_offset0; const double Ac_offset0 = fix->history_many[i][fix->c2r[ic]][18]; ddelta_bar[i] += wij*Ac_offset0/Acon0[i]*ddelta; // Multiply by 0.5 since displacement is shared equally between deformable particles. //std::cout << delta << ", " << delta_offset0 << " || " << Ac_offset0 << ", " << Acon0[i] << ", " << ddelta << ", " << ddelta_bar[i] << std::endl; } } - } + } } } @@ -536,4 +536,4 @@ void FixMDRmeanSurfDisp::pre_force(int) //std::cout << radius[i] << ", " << dR << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << dRdenominator[i] - 4.0*M_PI*pow(R,2.0) << std::endl; //std::cout << "Fix radius update setup has been entered !!!" << std::endl; -//std::cout << Ro[0] << ", " << Vgeo[0] << ", " << Velas[0] << std::endl; \ No newline at end of file +//std::cout << Ro[0] << ", " << Vgeo[0] << ", " << Velas[0] << std::endl; diff --git a/src/fix_mdr_mean_surf_disp.h b/src/fix_mdr_mean_surf_disp.h index be5be15907..daacb660e0 100644 --- a/src/fix_mdr_mean_surf_disp.h +++ b/src/fix_mdr_mean_surf_disp.h @@ -25,10 +25,10 @@ FixStyle(mdr/mean/surf/disp,FixMDRmeanSurfDisp); namespace LAMMPS_NS { class FixMDRmeanSurfDisp : public Fix { - public: - double * Acon0; + public: + double * Acon0; double * ddelta_bar; - + FixMDRmeanSurfDisp(class LAMMPS *, int, char **); int setmask() override; void setup_pre_force(int) override; diff --git a/src/fix_mdr_radius_update.cpp b/src/fix_mdr_radius_update.cpp index 010fad64d5..3f226fb87a 100644 --- a/src/fix_mdr_radius_update.cpp +++ b/src/fix_mdr_radius_update.cpp @@ -35,7 +35,7 @@ #include "gran_sub_mod_normal.h" #include "update.h" #include "comm.h" -#include +#include #include using namespace LAMMPS_NS; @@ -69,19 +69,19 @@ void FixMDRradiusUpdate::setup_pre_force(int /*vflag*/) 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_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_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); + 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]; @@ -89,9 +89,9 @@ void FixMDRradiusUpdate::setup_pre_force(int /*vflag*/) eps_bar = atom->dvector[index_eps_bar]; dRnumerator = atom->dvector[index_dRnumerator]; dRdenominator = atom->dvector[index_dRdenominator]; - Acon0 = atom->dvector[index_Acon0]; + Acon0 = atom->dvector[index_Acon0]; Acon1 = atom->dvector[index_Acon1]; - Atot = atom->dvector[index_Atot]; + Atot = atom->dvector[index_Atot]; Atot_sum = atom->dvector[index_Atot_sum]; ddelta_bar = atom->dvector[index_ddelta_bar]; psi = atom->dvector[index_psi]; @@ -116,19 +116,19 @@ void FixMDRradiusUpdate::setup(int /*vflag*/) 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_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_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); + 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]; @@ -136,9 +136,9 @@ void FixMDRradiusUpdate::setup(int /*vflag*/) eps_bar = atom->dvector[index_eps_bar]; dRnumerator = atom->dvector[index_dRnumerator]; dRdenominator = atom->dvector[index_dRdenominator]; - Acon0 = atom->dvector[index_Acon0]; + Acon0 = atom->dvector[index_Acon0]; Acon1 = atom->dvector[index_Acon1]; - Atot = atom->dvector[index_Atot]; + Atot = atom->dvector[index_Atot]; Atot_sum = atom->dvector[index_Atot_sum]; ddelta_bar = atom->dvector[index_ddelta_bar]; psi = atom->dvector[index_psi]; @@ -170,20 +170,20 @@ void FixMDRradiusUpdate::pre_force(int) //std::cout << "Preforce was called radius update" << std::endl; - // assign correct value to initially non-zero MDR particle history variables + // assign correct value to initially non-zero MDR particle history variables //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_Atot = atom->find_custom("Atot",tmp1,tmp2); + //int index_Atot = atom->find_custom("Atot",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); + //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); //double * Ro = atom->dvector[index_Ro]; //double * Vgeo = atom->dvector[index_Vgeo]; //double * Velas = atom->dvector[index_Velas]; @@ -196,13 +196,13 @@ void FixMDRradiusUpdate::pre_force(int) //double * history_setup_flag = atom->dvector[index_history_setup_flag]; //double * contacts = atom->dvector[index_contacts]; //double * adhesive_length = atom->dvector[index_adhesive_length]; - + double *radius = atom->radius; int nlocal = atom->nlocal; //std::cout << "ntotal " << ntotal << ", nlocal " << atom->nlocal << ", nghost " << atom->nghost << std::endl; - for (int i = 0; i < nlocal; i++) { + for (int i = 0; i < nlocal; i++) { if (history_setup_flag[i] < 1e-16) { Ro[i] = radius[i]; Vgeo[i] = 4.0/3.0*M_PI*pow(Ro[i],3.0); @@ -232,7 +232,7 @@ int FixMDRradiusUpdate::pack_forward_comm(int n, int *list, double *buf, int /*p buf[m++] = Vgeo[j]; // 2 buf[m++] = Velas[j]; // 3 buf[m++] = Vcaps[j]; // 4 - buf[m++] = eps_bar[j]; // 5 + buf[m++] = eps_bar[j]; // 5 buf[m++] = dRnumerator[j]; // 6 buf[m++] = dRdenominator[j]; // 7 buf[m++] = Acon0[j]; // 8 @@ -261,7 +261,7 @@ void FixMDRradiusUpdate::unpack_forward_comm(int n, int first, double *buf) Vgeo[i] = buf[m++]; // 2 Velas[i] = buf[m++]; // 3 Vcaps[i] = buf[m++]; // 4 - eps_bar[i] = buf[m++]; // 5 + eps_bar[i] = buf[m++]; // 5 dRnumerator[i] = buf[m++]; // 6 dRdenominator[i] = buf[m++]; // 7 Acon0[i] = buf[m++]; // 8 @@ -282,7 +282,7 @@ void FixMDRradiusUpdate::unpack_forward_comm(int n, int first, double *buf) void FixMDRradiusUpdate::end_of_step() { - + //comm->forward_comm(this); //std::cout << "end_of_step() was called radius update" << std::endl; @@ -293,13 +293,13 @@ void FixMDRradiusUpdate::end_of_step() int nlocal = atom->nlocal; double sigmaxx_sum = 0.0; double sigmayy_sum = 0.0; - double sigmazz_sum = 0.0; + double sigmazz_sum = 0.0; double Vparticles = 0.0; //std::cout << "New step" << std::endl; for (int i = 0; i < nlocal; i++) { - + const double R = radius[i]; Atot[i] = 4.0*M_PI*pow(R,2.0) + Atot_sum[i]; @@ -316,7 +316,7 @@ void FixMDRradiusUpdate::end_of_step() //} const double dR = std::max(dRnumerator[i]/(dRdenominator[i] - 4.0*M_PI*pow(R,2.0)),0.0); - if (psi_b[i] < psi[i]) { + if (psi_b[i] < psi[i]) { if ((radius[i] + dR) < (1.5*Ro[i])) radius[i] += dR; //radius[i] += dR; } @@ -327,10 +327,10 @@ void FixMDRradiusUpdate::end_of_step() //if (atom->tag[i] == 9){ // std::cout << i << ", radius: " << radius[i] << ", dR: " << dR << ", dRnum: " << dRnumerator[i] << ", dRdenom " << dRdenominator[i] << ", dRdem_full " << dRdenominator[i] - 4.0*M_PI*pow(R,2.0) << std::endl; //} - - + + //std::cout << i << ", " << radius[i] << " | " << psi_b[i] << ", " << psi[i] << " | " << Acon1[i] << ", " << Atot[i] << std::endl; - + Velas[i] = Vo*(1.0 + eps_bar[i]); Vcaps[i] = 0.0; @@ -339,7 +339,7 @@ void FixMDRradiusUpdate::end_of_step() dRdenominator[i] = 0.0; Acon0[i] = Acon1[i]; Acon1[i] = 0.0; - //std::cout << "Acon reset: " << Acon0[i] << ", " << Acon1[i] << std::endl; + //std::cout << "Acon reset: " << Acon0[i] << ", " << Acon1[i] << std::endl; Atot_sum[i] = 0.0; ddelta_bar[i] = 0.0; //adhesive_length[i] = adhesive_length[i]/contacts[i]; // convert adhesive length to average aAdh for each particle @@ -381,4 +381,4 @@ void FixMDRradiusUpdate::end_of_step() //std::cout << radius[i] << ", " << dR << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << dRdenominator[i] - 4.0*M_PI*pow(R,2.0) << std::endl; //std::cout << "Fix radius update setup has been entered !!!" << std::endl; -//std::cout << Ro[0] << ", " << Vgeo[0] << ", " << Velas[0] << std::endl; \ No newline at end of file +//std::cout << Ro[0] << ", " << Vgeo[0] << ", " << Velas[0] << std::endl; diff --git a/src/fix_mdr_radius_update.h b/src/fix_mdr_radius_update.h index e53366c25e..f431c0d13d 100644 --- a/src/fix_mdr_radius_update.h +++ b/src/fix_mdr_radius_update.h @@ -31,7 +31,7 @@ class FixMDRradiusUpdate : public Fix { double * Velas; double * Vcaps; double * eps_bar; - double * dRnumerator; + double * dRnumerator; double * dRdenominator; double * Acon0; double * Acon1; @@ -43,10 +43,10 @@ class FixMDRradiusUpdate : public Fix { double * sigmaxx; double * sigmayy; double * sigmazz; - double * history_setup_flag; + double * history_setup_flag; double * contacts; double * adhesive_length; - + FixMDRradiusUpdate(class LAMMPS *, int, char **); int setmask() override; void setup(int) override; diff --git a/src/variable.cpp b/src/variable.cpp index 350d6a7141..b581e4ce0d 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -3661,7 +3661,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree, Tree **tree strcmp(word,"logfreq") != 0 && strcmp(word,"logfreq2") != 0 && strcmp(word,"logfreq3") != 0 && strcmp(word,"stride") != 0 && strcmp(word,"stride2") != 0 && strcmp(word,"vdisplace") != 0 && - strcmp(word,"swiggle") != 0 && strcmp(word,"cwiggle") != 0 && strcmp(word,"sign") != 0) + strcmp(word,"swiggle") != 0 && strcmp(word,"cwiggle") != 0 && strcmp(word,"sign") != 0) return 0; // parse contents for comma-separated args