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     (if not contributing author is listed, this file has been contributed
36     by the core developer)
37 
38     Copyright 2012-     DCS Computing GmbH, Linz
39     Copyright 2009-2012 JKU Linz
40     Copyright 2016-     CFDEMresearch GmbH, Linz
41 ------------------------------------------------------------------------- */
42 
43 #ifdef PAIR_CLASS
44 
45 #else
46 
47 #ifndef LMP_PAIR_GRAN_H
48 #define LMP_PAIR_GRAN_H
49 
50 #include "pair.h"
51 #include "atom.h"
52 #include "region.h"
53 #include "compute_pair_gran_local.h"
54 #include "contact_interface.h"
55 #include "fix_relax_contacts.h"
56 #include "fix_property_atom.h"
57 #include <vector>
58 #include <string>
59 
60 namespace LCM = LIGGGHTS::ContactModels;
61 
62 namespace LAMMPS_NS {
63 
64 class ComputePairGranLocal;
65 
66 class PairGran : public Pair, public LIGGGHTS::IContactHistorySetup {
67 public:
68 
69   friend class FixWallGran;
70   friend class FixCheckTimestepGran;
71 
72   PairGran(class LAMMPS *);
73   ~PairGran();
74 
75   /* INHERITED FROM Pair */
76 
77   virtual void compute(int eflag, int vflag);
78   virtual void compute_pgl(int eflag, int vflag);
79   virtual void settings(int, char **) = 0;
80   virtual void coeff(int, char **);
81   virtual void init_style();
init_granular()82   virtual void init_granular() {}
83   virtual void init_list(int, class NeighList *);
84   virtual double init_one(int, int);
85   int pack_comm(int n, int *list,double *buf, int pbc_flag, int *pbc);
86   void unpack_comm(int n, int first, double *buf);
87   virtual void write_restart(FILE *);
88   virtual void read_restart(FILE *, const int major, const int minor);
write_restart_settings(FILE *)89   virtual void write_restart_settings(FILE *){}
read_restart_settings(FILE *,const int major,const int minor)90   virtual void read_restart_settings(FILE *, const int major, const int minor){}
91   virtual bool contact_match(const std::string mtype, const std::string model) = 0;
92   virtual void reset_dt();
93   double memory_usage();
94 
95   virtual int64_t hashcode() = 0;
96 
cplenable()97   int  cplenable()
98   { return cpl_enable; }
99 
100   void register_compute_pair_local(class ComputePairGranLocal *,int&);
101   void unregister_compute_pair_local(class ComputePairGranLocal *ptr);
102 
103   void cpl_add_pair(LCM::SurfacesIntersectData & sidata, LCM::ForceData & i_forces);
104 
105   void cpl_pair_finalize();
106 
107   /* PUBLIC ACCESS FUNCTIONS */
108 
is_history()109   int is_history()
110   { return history; }
111 
dnum_pair()112   int dnum_pair()
113   { return dnum_pairgran; }
114 
dnum()115   inline int dnum()
116   { return dnum_all; }
117 
cpl()118   inline class ComputePairGranLocal * cpl() const
119   { return cpl_; }
120 
storeContactForces()121   inline bool storeContactForces() const
122   { return store_contact_forces_; }
123 
storeContactForcesEvery()124   inline int storeContactForcesEvery() const
125   { return store_contact_forces_every_; }
126 
storeContactForcesStress()127   inline bool storeContactForcesStress()
128   { return store_contact_forces_stress_; }
129 
storeSumDelta()130   inline bool storeSumDelta()
131   { return store_multicontact_data_; }
132 
freeze_group_bit()133   inline int freeze_group_bit() const
134   { return freeze_group_bit_; }
135 
computeflag()136   inline int computeflag() const
137   { return computeflag_; }
138 
shearupdate()139   inline int shearupdate() const
140   { return shearupdate_; }
141 
fix_contact_forces()142   class FixContactPropertyAtom * fix_contact_forces() const
143   { return fix_contact_forces_;  }
144 
fix_contact_forces_stress()145   class FixContactPropertyAtom * fix_contact_forces_stress()
146   { return fix_contact_forces_stress_; }
147 
fix_store_multicontact_delta()148   class FixContactPropertyAtom * fix_store_multicontact_delta()
149   { return fix_store_multicontact_data_; }
150 
fr_pair()151   class FixRigid* fr_pair() const
152   { return fix_rigid; }
153 
mr_pair()154   double * mr_pair() const
155   { return mass_rigid; }
156 
relax(int i)157   double relax(int i)
158   { return (fix_relax_ ? fix_relax_->factor_relax(i) : 1.); }
159 
160   virtual double stressStrainExponent() = 0;
161 
162   int fix_extra_dnum_index(class Fix *fix);
163 
164   void *extract(const char *str, int &dim);
165 
add_history_value(std::string name,std::string newtonflag)166   int add_history_value(std::string name, std::string newtonflag)
167   {
168     int offset = history_arg.size();
169     history = true;
170     history_arg.push_back(HistoryArg(name, newtonflag));
171     dnum_pairgran++;
172     return offset;
173   }
174 
add_dissipated_energy(const double e)175   void add_dissipated_energy(const double e)
176   { dissipated_energy_ += e; }
177 
get_dissipated_energy()178   double get_dissipated_energy()
179   { return dissipated_energy_; }
180 
do_store_contact_forces()181   void do_store_contact_forces()
182   { store_contact_forces_ = true; }
183 
do_store_contact_forces_every(int ev)184   void do_store_contact_forces_every(int ev)
185   { store_contact_forces_ = true;  store_contact_forces_every_ = ev; }
186 
do_relax_region(FixRelaxContacts * _fr)187   void do_relax_region(FixRelaxContacts *_fr)
188   { fix_relax_ = _fr; }
189 
do_store_contact_forces_stress()190   void do_store_contact_forces_stress()
191   { store_contact_forces_stress_ = true; }
192 
do_store_multicontact_data()193   void do_store_multicontact_data()
194   { store_multicontact_data_ = true; }
195 
get_fix_history()196   class FixContactHistory* get_fix_history() const
197   { return fix_history; }
198 
store_sum_normal_force()199   bool store_sum_normal_force() const
200   { return fix_sum_normal_force_ != NULL; }
201 
get_sum_normal_force_ptr(const int i)202   double * get_sum_normal_force_ptr(const int i)
203   { return &(fix_sum_normal_force_->vector_atom[i]); }
204 
205  protected:
206 
207   struct HistoryArg {
208     std::string name;
209     std::string newtonflag;
210 
HistoryArgHistoryArg211     HistoryArg(std::string name, std::string newtonflag) : name(name), newtonflag(newtonflag) {}
212   };
213 
214   std::vector<HistoryArg> history_arg;
215 
history_args(char ** args)216   void history_args(char ** args) {
217     for(size_t i = 0; i < history_arg.size(); i++) {
218       args[2*i] = (char*)history_arg[i].name.c_str();
219       args[2*i+1] = (char*)history_arg[i].newtonflag.c_str();
220     }
221   }
222 
223   virtual void compute_force(int eflag, int vflag,int addflag) = 0;
224   virtual bool forceoff();
225 
226   virtual void updatePtrs();
227 
228   // for parsing settings() args
229   int iarg_;
230 
231   char * suffix;
232   int neighprev;
233 
234   // stuff for tracking energy
235   int energytrack_enable;
236   class FixPropertyAtom* fppaCPEn; //collision potential energy normal
237   class FixPropertyAtom* fppaCDEn; //collision dissipation energy normal
238   class FixPropertyAtom* fppaCPEt; //collision potential energy tang
239   class FixPropertyAtom* fppaCDEVt; //collision dissipation energy viscous tang
240   class FixPropertyAtom* fppaCDEFt; //collision dissipation energy friction tang
241   class FixPropertyAtom* fppaCTFW; //collision tangential force work
242   class FixPropertyAtom* fppaDEH; //dissipation energy of history term (viscous and friction, accumulated over time)
243   double *CPEn, *CDEn, *CPEt, *CDEVt, *CDEFt, *CTFW, *DEH;
244 
245   // stuff for compute pair gran local
246   int cpl_enable;
247   class ComputePairGranLocal *cpl_;
248 
249   // storage for per-contact forces and torque
250   bool store_contact_forces_;
251   int store_contact_forces_every_;
252   class FixContactPropertyAtom *fix_contact_forces_;
253 
254   // storage for per-contact forces and relative position (for goldhirsch stress model)
255   bool store_contact_forces_stress_;
256   class FixContactPropertyAtom *fix_contact_forces_stress_;
257   // storage for per contact delta data (for multicontact models)
258   bool store_multicontact_data_;
259   class FixContactPropertyAtom *fix_store_multicontact_data_;
260   // storage for simplistic pressure computation via normal forces
261   class FixPropertyAtom *fix_sum_normal_force_;
262 
263   // storage of rigid body masses for use in granular interactions
264 
265   class FixRigid *fix_rigid; // ptr to rigid body fix, NULL if none
266   double *mass_rigid;        // rigid mass for owned+ghost atoms
267   int nmax;                  // allocated size of mass_rigid
268 
269   double dt;
270   int freeze_group_bit_;
271 
272   // contact history
273   int history;
274   int dnum_pairgran;
275   class FixContactHistory *fix_history;
276   int shearupdate_;
277 
278   int computeflag_;
279 
280   double *onerad_dynamic,*onerad_frozen;
281   double *maxrad_dynamic,*maxrad_frozen;
282 
283   bool needs_neighlist;
284 
285   void allocate();
286 
287  private:
288 
289   // shear history
290   int dnum_all;
291 
292   int nfix;
293   Fix **fix_dnum;
294   int *dnum_index;
295 
296   FixRelaxContacts *fix_relax_;
297 
298   // dissipated energy in wall -> particle contacts
299   double dissipated_energy_;
300 };
301 
302 }
303 
304 #endif
305 #endif
306