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 Arno Mayrhofer (CFDEMresearch GmbH, Linz) 40 41 Copyright 2012- DCS Computing GmbH, Linz 42 Copyright 2009-2012 JKU Linz 43 Copyright 2016- CFDEMresearch GmbH, Linz 44 ------------------------------------------------------------------------- */ 45 46 #ifndef LMP_FIX_WALL_GRAN_BASE_H 47 #define LMP_FIX_WALL_GRAN_BASE_H 48 49 #include "fix_wall_gran.h" 50 #include "fix_contact_property_atom_wall.h" 51 #include "contact_interface.h" 52 #include "compute_pair_gran_local.h" 53 #include "tri_mesh.h" 54 #include "settings.h" 55 #include <string.h> 56 #include "force.h" 57 #include <stdlib.h> 58 #include "contact_models.h" 59 #include "granular_wall.h" 60 #include "fix_calculate_energy_wall.h" 61 62 #ifdef SUPERQUADRIC_ACTIVE_FLAG 63 #include "math_const.h" 64 #endif 65 66 namespace LIGGGHTS { 67 using namespace ContactModels; 68 namespace Walls { 69 70 template<typename ContactModel> 71 class Granular : private Pointers, public IGranularWall { 72 ContactModel cmodel; 73 FixWallGran * parent; 74 int dissipation_offset_; 75 FixCalculateWallEnergy *fix_dissipated_energy_; 76 77 public: Granular(LAMMPS * lmp,FixWallGran * parent,const int64_t hash)78 Granular(LAMMPS * lmp, FixWallGran * parent, const int64_t hash) : 79 Pointers(lmp), 80 cmodel(lmp, parent,true /*is_wall*/, hash), 81 parent(parent), 82 dissipation_offset_(-1), 83 fix_dissipated_energy_(NULL) 84 { 85 } 86 ~Granular()87 virtual ~Granular() 88 { 89 } 90 init_granular()91 virtual void init_granular() { 92 cmodel.connectToProperties(force->registry); 93 94 #ifdef LIGGGHTS_DEBUG 95 if(comm->me == 0) { 96 fprintf(screen, "==== WALL %s GLOBAL PROPERTIES ====\n", parent->id); 97 force->registry.print_all(screen); 98 fprintf(screen, "==== WALL %s GLOBAL PROPERTIES ====\n", parent->id); 99 100 fprintf(logfile, "==== WALL %s GLOBAL PROPERTIES ====\n", parent->id); 101 force->registry.print_all(logfile); 102 fprintf(logfile, "==== WALL %s GLOBAL PROPERTIES ====\n", parent->id); 103 } 104 #endif 105 } 106 settings(int nargs,char ** args,IContactHistorySetup * hsetup)107 virtual void settings(int nargs, char ** args, IContactHistorySetup *hsetup) { 108 Settings settings(lmp); 109 cmodel.registerSettings(settings); 110 bool success = settings.parseArguments(nargs, args); 111 cmodel.postSettings(hsetup); 112 113 #ifdef LIGGGHTS_DEBUG 114 if(comm->me == 0) { 115 fprintf(screen, "==== WALL %s SETTINGS ====\n", parent->id); 116 settings.print_all(screen); 117 fprintf(screen, "==== WALL %s SETTINGS ====\n", parent->id); 118 119 fprintf(logfile, "==== WALL %s SETTINGS ====\n", parent->id); 120 settings.print_all(logfile); 121 fprintf(logfile, "==== WALL %s SETTINGS ====\n", parent->id); 122 } 123 #endif 124 125 dissipation_offset_ = get_history_offset("dissipation_force"); 126 fix_dissipated_energy_ = static_cast<FixCalculateWallEnergy*>(modify->find_fix_style("calculate/wall_dissipated_energy", 0)); 127 if (dissipation_offset_ >= 0 && !fix_dissipated_energy_) 128 error->one(FLERR, "Could not find fix calculate/wall_dissipated_energy"); 129 130 if(!success) { 131 error->fix_error(FLERR, parent, settings.error_message.c_str()); 132 } 133 } 134 force_update(double * const f,double * const torque,const ForceData & forces)135 inline void force_update(double * const f, double * const torque, 136 const ForceData & forces) { 137 for (int coord = 0; coord < 3; coord++) { 138 f[coord] += forces.delta_F[coord]; 139 torque[coord] += forces.delta_torque[coord]; 140 } 141 } 142 checkSurfaceIntersect(SurfacesIntersectData & sidata)143 virtual bool checkSurfaceIntersect(SurfacesIntersectData & sidata) 144 { 145 return cmodel.checkSurfaceIntersect(sidata); 146 } 147 148 virtual void compute_force(FixWallGran * wg, SurfacesIntersectData & sidata, bool intersectflag,double *vwall, class FixMeshSurface * fix_mesh = 0, int iMesh = 0, class TriMesh *mesh = 0,int iTri = 0) 149 { 150 const int ip = sidata.i; 151 152 double *x = atom->x[ip]; 153 double *f = atom->f[ip]; 154 double *torque = atom->torque[ip]; 155 double *v = atom->v[ip]; 156 double *omega = atom->omega[ip]; 157 double mass = atom->rmass[ip]; 158 int *type = atom->type; 159 160 #ifdef LIGGGHTS_DEBUG 161 if(std::isnan(vectorMag3D(x))) 162 error->one(FLERR,"x is NaN!"); 163 if(std::isnan(vectorMag3D(f))) 164 error->one(FLERR,"f is NaN!"); 165 if(std::isnan(vectorMag3D(torque))) 166 error->one(FLERR,"torque is NaN!"); 167 if(std::isnan(vectorMag3D(omega))) 168 error->one(FLERR,"omega is NaN!"); 169 #endif 170 171 // copy collision data to struct (compiler can figure out a better way to 172 // interleave these stores with the double calculations above. 173 ForceData i_forces; 174 ForceData j_forces; 175 memset(&i_forces, 0, sizeof(ForceData)); 176 memset(&j_forces, 0, sizeof(ForceData)); 177 sidata.v_i = v; 178 sidata.v_j = vwall; 179 sidata.omega_i = omega; 180 181 sidata.r = sidata.radi - sidata.deltan; // sign corrected, because negative value is passed 182 sidata.rsq = sidata.r*sidata.r; 183 const double rinv = 1.0/sidata.r; 184 sidata.rinv = rinv; 185 sidata.area_ratio = 1.; 186 187 //store type here as negative in case of primitive wall! 188 sidata.j = mesh ? iTri : -wg->atom_type_wall(); 189 sidata.contact_flags = NULL; 190 sidata.itype = type[ip]; 191 192 if(wg->fix_rigid() && wg->body(ip) >= 0) 193 mass = wg->masstotal(wg->body(ip)); 194 195 sidata.meff = mass; 196 sidata.mi = mass; 197 198 sidata.computeflag = wg->computeflag(); 199 sidata.shearupdate = wg->shearupdate(); 200 sidata.jtype = wg->atom_type_wall(); 201 202 double force_old[3]={}, f_pw[3]; 203 204 // if force should be stored - remember old force 205 if(wg->store_force() || fix_mesh) 206 vectorCopy3D(f,force_old); 207 208 // add to cwl 209 if(wg->compute_wall_gran_local() && wg->addflag()) 210 { 211 double contactPoint[3]; 212 vectorSubtract3D(x,sidata.delta,contactPoint); 213 #ifdef SUPERQUADRIC_ACTIVE_FLAG 214 if(atom->superquadric_flag) { 215 vectorCopy3D(sidata.contact_point, contactPoint); 216 } 217 #endif 218 wg->compute_wall_gran_local()->add_wall_1(iMesh,mesh->id(iTri),ip,contactPoint,vwall); 219 } 220 221 #ifdef SUPERQUADRIC_ACTIVE_FLAG 222 double enx, eny, enz; 223 if(atom->superquadric_flag) { 224 #ifdef LIGGGHTS_DEBUG 225 if(vectorMag3D(sidata.delta) == 0.0) 226 error->one(FLERR, "vectorMag3D(sidata.delta) == 0.0"); 227 #endif 228 const double delta_inv = 1.0 / vectorMag3D(sidata.delta); 229 #ifdef LIGGGHTS_DEBUG 230 if(std::isnan(delta_inv)) 231 error->one(FLERR, "delta_inv is NaN!"); 232 #endif 233 enx = sidata.delta[0] * delta_inv; 234 eny = sidata.delta[1] * delta_inv; 235 enz = sidata.delta[2] * delta_inv; 236 sidata.radi = cbrt(0.75 * atom->volume[ip] / M_PI); 237 #ifdef LIGGGHTS_DEBUG 238 if(std::isnan(sidata.radi)) 239 error->one(FLERR, "delta_inv is NaN!"); 240 #endif 241 } else { // sphere case 242 enx = sidata.delta[0] * rinv; 243 eny = sidata.delta[1] * rinv; 244 enz = sidata.delta[2] * rinv; 245 } 246 #else // sphere case 247 const double enx = sidata.delta[0] * rinv; 248 const double eny = sidata.delta[1] * rinv; 249 const double enz = sidata.delta[2] * rinv; 250 #endif 251 252 sidata.radsum = sidata.radi; 253 254 #ifdef CONVEX_ACTIVE_FLAG 255 if (!atom->shapetype_flag) 256 #endif 257 { 258 sidata.en[0] = enx; 259 sidata.en[1] = eny; 260 sidata.en[2] = enz; 261 } 262 263 double delta[3]; 264 if (dissipation_offset_ >= 0 && sidata.computeflag && sidata.shearupdate) 265 { 266 sidata.fix_mesh->triMesh()->get_global_vel(delta); 267 vectorScalarMult3D(delta, update->dt); 268 // displacement force from the previous time step 269 double * const diss_force = &sidata.contact_history[dissipation_offset_]; 270 // in the scalar product with the current mesh delta 271 const double dissipated_energy = vectorDot3D(delta, diss_force)*0.5; 272 fix_dissipated_energy_->dissipate_energy(dissipated_energy, true); 273 // reset the dissipation force 274 diss_force[0] = 0.0; 275 diss_force[1] = 0.0; 276 diss_force[2] = 0.0; 277 } 278 279 if (intersectflag) 280 { 281 cmodel.surfacesIntersect(sidata, i_forces, j_forces); 282 cmodel.endSurfacesIntersect(sidata, mesh, i_forces, j_forces); 283 // if there is a surface touch, there will always be a force 284 sidata.has_force_update = true; 285 } 286 // surfacesClose is not supported for convex particles 287 else if (!atom->shapetype_flag) 288 { 289 // apply force update only if selected contact models have requested it 290 sidata.has_force_update = false; 291 cmodel.surfacesClose(sidata, i_forces, j_forces); 292 } 293 294 if (sidata.is_wall && dissipation_offset_ >= 0 && sidata.computeflag && sidata.shearupdate) 295 { 296 // new displacement force from this time step 297 double * const diss_force = &sidata.contact_history[dissipation_offset_]; 298 // in the scalar product with the mesh displacement from earlier on 299 const double dissipated_energy = vectorDot3D(delta, diss_force)*0.5; 300 //printf("pdis %e %e\n", update->get_cur_time(), dissipated_energy); 301 fix_dissipated_energy_->dissipate_energy(dissipated_energy, false); 302 } 303 304 if(sidata.computeflag) 305 { 306 if (sidata.has_force_update) 307 force_update(f, torque, i_forces); 308 // summation of f.n to compute a simplistic pressure 309 if (wg->store_sum_normal_force()) 310 { 311 double * const iforce = wg->get_sum_normal_force_ptr(ip); 312 const double fDotN = vectorDot3D(i_forces.delta_F, sidata.en); 313 *iforce += fDotN; 314 } 315 } 316 317 if (wg->store_force_contact() && 0 == update->ntimestep % wg->store_force_contact_every()) 318 { 319 wg->add_contactforce_wall(ip,i_forces,mesh?mesh->id(iTri):0); 320 } 321 322 if (wg->store_force_contact_stress()) 323 wg->add_contactforce_stress_wall(ip, i_forces, sidata.delta, vwall, mesh?mesh->id(iTri):0); 324 325 if(wg->compute_wall_gran_local() && wg->addflag()) 326 { 327 const double fx = i_forces.delta_F[0]; 328 const double fy = i_forces.delta_F[1]; 329 const double fz = i_forces.delta_F[2]; 330 const double tor1 = i_forces.delta_torque[0]*sidata.area_ratio; 331 const double tor2 = i_forces.delta_torque[1]*sidata.area_ratio; 332 const double tor3 = i_forces.delta_torque[2]*sidata.area_ratio; 333 double normal[3]; 334 vectorCopy3D(sidata.en, normal); 335 vectorNegate3D(normal); 336 wg->compute_wall_gran_local()->add_wall_2(sidata.i,fx,fy,fz,tor1,tor2,tor3,sidata.contact_history,sidata.rsq, normal); 337 } 338 339 // add heat flux 340 341 if(wg->heattransfer_flag()) 342 wg->addHeatFlux(mesh,ip,sidata.radi,sidata.deltan,1.); 343 344 // if force should be stored or evaluated 345 if(sidata.has_force_update && (wg->store_force() || fix_mesh) ) 346 { 347 vectorSubtract3D(f,force_old,f_pw); 348 349 if(wg->store_force()) 350 vectorAdd3D (wg->fix_wallforce()->array_atom[ip], f_pw, wg->fix_wallforce()->array_atom[ip]); 351 352 if (fix_mesh) 353 { 354 double delta[3]; 355 delta[0] = -sidata.delta[0]; 356 delta[1] = -sidata.delta[1]; 357 delta[2] = -sidata.delta[2]; 358 fix_mesh->add_particle_contribution (ip,f_pw,delta,iTri,vwall); 359 } 360 } 361 } 362 get_history_offset(const std::string hname)363 int get_history_offset(const std::string hname) 364 { 365 return cmodel.get_history_offset(hname); 366 } 367 contact_match(const std::string mtype,const std::string model)368 bool contact_match(const std::string mtype, const std::string model) 369 { 370 return cmodel.contact_match(mtype, model); 371 } 372 373 }; 374 375 } 376 377 } 378 379 #endif 380