1 /* ---------------------------------------------------------------------- 2 This is the 3 4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗ 5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝ 6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗ 7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║ 8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║ 9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝® 10 11 DEM simulation engine, released by 12 DCS Computing Gmbh, Linz, Austria 13 http://www.dcs-computing.com, office@dcs-computing.com 14 15 LIGGGHTS® is part of CFDEM®project: 16 http://www.liggghts.com | http://www.cfdem.com 17 18 Core developer and main author: 19 Christoph Kloss, christoph.kloss@dcs-computing.com 20 21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public 22 License, version 2 or later. It is distributed in the hope that it will 23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have 25 received a copy of the GNU General Public License along with LIGGGHTS®. 26 If not, see http://www.gnu.org/licenses . See also top-level README 27 and LICENSE files. 28 29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH, 30 the producer of the LIGGGHTS® software and the CFDEM®coupling software 31 See http://www.cfdem.com/terms-trademark-policy for details. 32 33 ------------------------------------------------------------------------- 34 Contributing author and copyright for this file: 35 36 Christoph Kloss (DCS Computing GmbH, Linz) 37 Christoph Kloss (JKU Linz) 38 Richard Berger (JKU Linz) 39 40 Copyright 2012- DCS Computing GmbH, Linz 41 Copyright 2009-2012 JKU Linz 42 ------------------------------------------------------------------------- */ 43 44 #ifdef SURFACE_MODEL 45 SURFACE_MODEL(SURFACE_DEFAULT,default,0) 46 #else 47 #ifndef SURFACE_MODEL_DEFAULT_H_ 48 #define SURFACE_MODEL_DEFAULT_H_ 49 #include "contact_models.h" 50 #include <cmath> 51 #include "atom.h" 52 #include "force.h" 53 #include "update.h" 54 #include "modify.h" 55 #include "fix_property_atom.h" 56 #include "surface_model_base.h" 57 58 namespace LIGGGHTS { 59 namespace ContactModels 60 { 61 template<> 62 class SurfaceModel<SURFACE_DEFAULT> : public SurfaceModelBase 63 { 64 public: 65 SurfaceModel(LAMMPS * lmp, IContactHistorySetup * hsetup, class ContactModelBase * cmb) : 66 SurfaceModelBase(lmp, hsetup, cmb), 67 elasticpotflag_(false), 68 dissipatedflag_(false), 69 delta_offset_(-1), 70 dissipation_offset_(-1), 71 fix_dissipated_(NULL) 72 { 73 74 } 75 76 inline void registerSettings(Settings& settings) 77 { 78 settings.registerOnOff("computeElasticPotential", elasticpotflag_, false); 79 settings.registerOnOff("computeDissipatedEnergy", dissipatedflag_, false); 80 } 81 82 inline void postSettings(IContactHistorySetup * hsetup, ContactModelBase *cmb) 83 { 84 if (dissipatedflag_) 85 { 86 if (cmb->is_wall()) 87 { 88 fix_dissipated_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("dissipated_energy_wall", "property/atom", "vector", 0, 0, "dissipated energy")); 89 if (!fix_dissipated_) 90 error->one(FLERR, "Could not find dissipated_energy_wall atom property. Ensure that fix calculate/wall_dissipated_energy is before fix wall/gran"); 91 } 92 else 93 { 94 char * fixarg[15]; 95 fixarg[0] = (char*)"dissipated_energy_"; 96 fixarg[1] = (char*)"all"; 97 fixarg[2] = (char*)"property/atom"; 98 fixarg[3] = (char*)"dissipated_energy"; 99 fixarg[4] = (char*)"vector"; 100 fixarg[5] = (char*)"yes"; 101 fixarg[6] = (char*)"yes"; 102 fixarg[7] = (char*)"no"; 103 fixarg[8] = (char*)"0.0"; // energy 104 fixarg[9] = (char*)"0.0"; // fx 105 fixarg[10] = (char*)"0.0"; // fy 106 fixarg[11] = (char*)"0.0"; // fz 107 fixarg[12] = (char*)"0.0"; // tx 108 fixarg[13] = (char*)"0.0"; // ty 109 fixarg[14] = (char*)"0.0"; // tz 110 fix_dissipated_ = modify->add_fix_property_atom(15, static_cast<char**>(fixarg), "dissipated energy"); 111 } 112 } 113 if (cmb->is_wall() && (dissipatedflag_ || elasticpotflag_)) 114 { 115 delta_offset_ = hsetup->add_history_value("delta_0", "1"); 116 hsetup->add_history_value("delta_1", "1"); 117 hsetup->add_history_value("delta_2", "1"); 118 cmb->add_history_offset("delta", delta_offset_); 119 if (dissipatedflag_) 120 { 121 dissipation_offset_ = hsetup->add_history_value("diss_f_0", "1"); 122 hsetup->add_history_value("diss_f_1", "1"); 123 hsetup->add_history_value("diss_f_2", "1"); 124 cmb->add_history_offset("dissipation_force", dissipation_offset_); 125 } 126 } 127 } 128 129 inline void connectToProperties(PropertyRegistry&) {} 130 131 inline bool checkSurfaceIntersect(SurfacesIntersectData & sidata) 132 { 133 #ifdef SUPERQUADRIC_ACTIVE_FLAG 134 sidata.is_non_spherical = false; 135 #endif 136 137 return true; 138 } 139 140 inline void surfacesIntersect(SurfacesIntersectData & sidata, ForceData&, ForceData&) 141 { 142 #ifdef SUPERQUADRIC_ACTIVE_FLAG 143 if(sidata.is_non_spherical) 144 error->one(FLERR,"Using default surface model for non-spherical particles!"); 145 #endif 146 const double enx = sidata.en[0]; 147 const double eny = sidata.en[1]; 148 const double enz = sidata.en[2]; 149 150 // relative translational velocity 151 const double vr1 = sidata.v_i[0] - sidata.v_j[0]; 152 const double vr2 = sidata.v_i[1] - sidata.v_j[1]; 153 const double vr3 = sidata.v_i[2] - sidata.v_j[2]; 154 155 // normal component 156 const double vn = vr1 * enx + vr2 * eny + vr3 * enz; 157 const double vn1 = vn * enx; 158 const double vn2 = vn * eny; 159 const double vn3 = vn * enz; 160 161 // tangential component 162 const double vt1 = vr1 - vn1; 163 const double vt2 = vr2 - vn2; 164 const double vt3 = vr3 - vn3; 165 166 // relative rotational velocity 167 const double deltan = sidata.radsum - sidata.r; 168 const double dx = sidata.delta[0]; 169 const double dy = sidata.delta[1]; 170 const double dz = sidata.delta[2]; 171 const double rinv = sidata.rinv; 172 double wr1, wr2, wr3; 173 174 if(sidata.is_wall) { 175 // in case of wall contact, r is the contact radius 176 const double cr = sidata.radi - 0.5*sidata.deltan; 177 wr1 = cr * sidata.omega_i[0] * rinv; 178 wr2 = cr * sidata.omega_i[1] * rinv; 179 wr3 = cr * sidata.omega_i[2] * rinv; 180 sidata.cri = cr; 181 } else { 182 const double cri = sidata.radi - 0.5 * deltan; 183 const double crj = sidata.radj - 0.5 * deltan; 184 wr1 = (cri * sidata.omega_i[0] + crj * sidata.omega_j[0]) * rinv; 185 wr2 = (cri * sidata.omega_i[1] + crj * sidata.omega_j[1]) * rinv; 186 wr3 = (cri * sidata.omega_i[2] + crj * sidata.omega_j[2]) * rinv; 187 sidata.cri = cri; 188 sidata.crj = crj; 189 } 190 191 // relative velocities 192 const double vtr1 = vt1 - (dz * wr2 - dy * wr3); 193 const double vtr2 = vt2 - (dx * wr3 - dz * wr1); 194 const double vtr3 = vt3 - (dy * wr1 - dx * wr2); 195 196 sidata.vn = vn; 197 sidata.deltan = deltan; 198 sidata.wr1 = wr1; 199 sidata.wr2 = wr2; 200 sidata.wr3 = wr3; 201 sidata.vtr1 = vtr1; 202 sidata.vtr2 = vtr2; 203 sidata.vtr3 = vtr3; 204 sidata.P_diss = 0.; 205 } 206 207 inline void endSurfacesIntersect(SurfacesIntersectData &sidata,TriMesh *, double * const) {} 208 inline void surfacesClose(SurfacesCloseData &scdata, ForceData&, ForceData&) {} 209 void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){} 210 void endPass(SurfacesIntersectData&, ForceData&, ForceData&){} 211 inline void tally_pp(double,int,int,int) {} 212 inline void tally_pw(double,int,int,int) {} 213 214 private: 215 bool elasticpotflag_; 216 bool dissipatedflag_; 217 int delta_offset_; 218 int dissipation_offset_; 219 FixPropertyAtom *fix_dissipated_; 220 }; 221 } 222 } 223 #endif // SURFACE_MODEL_DEFAULT_H_ 224 #endif 225