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