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     Alexander Podlozhnyuk (DCS Computing GmbH, Linz)
37 
38     Copyright 2014-     DCS Computing GmbH, Linz
39 ------------------------------------------------------------------------- */
40 
41 #ifdef SURFACE_MODEL
42 SURFACE_MODEL(SURFACE_SUPERQUADRIC,superquadric,5)
43 #else
44 #ifndef SURFACE_MODEL_SUPERQUADRIC_H_
45 #define SURFACE_MODEL_SUPERQUADRIC_H_
46 
47 #include "surface_model_base.h"
48 
49 #ifdef SUPERQUADRIC_ACTIVE_FLAG
50 
51 #include "contact_models.h"
52 #include <cmath>
53 #include <algorithm>
54 #include "atom.h"
55 #include "force.h"
56 #include "update.h"
57 #include "math_extra_liggghts_superquadric.h"
58 
59 namespace LIGGGHTS {
60 namespace ContactModels
61 {
62   template<>
63   class SurfaceModel<SURFACE_SUPERQUADRIC> : public SurfaceModelBase
64   {
65     int inequality_start_offset;
66     int particles_were_in_contact_offset;
67     int contact_point_offset;
68     int alpha1_offset;
69     int alpha2_offset;
70     int reff_offset;
71     Superquadric particle_i;
72     Superquadric particle_j;
73     enum {SURFACES_FAR, SURFACES_CLOSE, SURFACES_INTERSECT};
74   public:
75     SurfaceModel(LAMMPS * lmp, IContactHistorySetup* hsetup,class ContactModelBase *cmb) :
76         SurfaceModelBase(lmp, hsetup, cmb)
77     {
78       if(!atom->superquadric_flag)
79         error->one(FLERR,"Applying surface model superquadric to a non-superquadric particle!");
80       inequality_start_offset = hsetup->add_history_value("inequality_obb","0");
81       particles_were_in_contact_offset = hsetup->add_history_value("particles_were_in_contact","0");
82       contact_point_offset = hsetup->add_history_value("cpx", "0");
83       cmb->add_history_offset("contact_point_offset", contact_point_offset);
84       hsetup->add_history_value("cpy", "0");
85       hsetup->add_history_value("cpz", "0");
86       alpha1_offset = hsetup->add_history_value("a1", "0");
87       alpha2_offset = hsetup->add_history_value("a2", "0");
88       cmb->add_history_offset("alpha1_offset", alpha1_offset);
89       cmb->add_history_offset("alpha2_offset", alpha2_offset);
90       cmb->get_history_offset("alpha2_offset");
91     }
92 
93     inline void registerSettings(Settings& settings) {
94       settings.registerDoubleSetting("curvatureLimitFactor",curvatureLimitFactor, 0.0);
95       settings.registerYesNo("meanCurvature", meanCurvature, false);
96       settings.registerYesNo("gaussianCurvature", gaussianCurvature, false);
97       if(curvatureLimitFactor < 0.0)
98         error->one(FLERR,"Curvature limiter cannot be negative!");
99       if(!meanCurvature && !gaussianCurvature)
100         curvatureLimitFactor = 0.0;
101       if(meanCurvature && gaussianCurvature)
102         error->one(FLERR,"meanCurvature and gaussianCurvature cannot be simultaneously yes !");
103     }
104 
105     inline void postSettings(IContactHistorySetup * hsetup, ContactModelBase *cmb) {}
106     inline void connectToProperties(PropertyRegistry&) {}
107 
108     inline bool checkSurfaceIntersect(SurfacesIntersectData & sidata)
109     {
110       sidata.is_non_spherical = true;
111       bool particles_in_contact = false;
112       double *const prev_step_point = &sidata.contact_history[contact_point_offset]; //contact points
113       double *const inequality_start = &sidata.contact_history[inequality_start_offset];
114       double *const particles_were_in_contact = &sidata.contact_history[particles_were_in_contact_offset];
115 
116       const int iPart = sidata.i;
117       const int jPart = sidata.j;
118 
119 #ifdef LIGGGHTS_DEBUG
120       if(std::isnan(vectorMag3D(atom->x[iPart])))
121         error->one(FLERR,"atom->x[iPart] is NaN!");
122       if(std::isnan(vectorMag4D(atom->quaternion[iPart])))
123         error->one(FLERR,"atom->quaternion[iPart] is NaN!");
124       if(std::isnan(vectorMag3D(atom->x[jPart])))
125         error->one(FLERR,"atom->x[jPart] is NaN!");
126       if(std::isnan(vectorMag4D(atom->quaternion[jPart])))
127         error->one(FLERR,"atom->quaternion[jPart] is NaN!");
128 #endif
129       particle_i.set(atom->x[iPart], atom->quaternion[iPart], atom->shape[iPart], atom->blockiness[iPart]);
130       particle_j.set(atom->x[jPart], atom->quaternion[jPart], atom->shape[jPart], atom->blockiness[jPart]);
131 
132       unsigned int int_inequality_start = MathExtraLiggghtsNonspherical::round_int(*inequality_start);
133 
134       bool obb_intersect = false;
135       if(*particles_were_in_contact == SURFACES_INTERSECT)
136         obb_intersect = true; //particles had overlap on the previous time step, skipping OBB intersection check
137       else
138         obb_intersect = MathExtraLiggghtsNonspherical::obb_intersect(&particle_i, &particle_j, int_inequality_start);
139       if(obb_intersect) {//OBB intersect particles in possible contact
140 
141         double fi, fj;
142         const double ri = cbrt(particle_i.shape[0]*particle_i.shape[1]*particle_i.shape[2]);
143         const double rj = cbrt(particle_j.shape[0]*particle_j.shape[1]*particle_j.shape[2]);
144         double ratio = ri / (ri + rj);
145 
146         if(*particles_were_in_contact == SURFACES_FAR)
147           calc_contact_point_if_no_previous_point_avaialable(sidata, &particle_i, &particle_j, sidata.contact_point, fi, fj, this->error);
148         else
149           calc_contact_point_using_prev_step(sidata, &particle_i, &particle_j, ratio, update->dt, prev_step_point, sidata.contact_point, fi, fj, this->error);
150         vectorCopy3D(sidata.contact_point, prev_step_point); //store contact point in contact history for the next DEM time step
151 
152 #ifdef LIGGGHTS_DEBUG
153         if(std::isnan(vectorMag3D(sidata.contact_point)))
154           error->one(FLERR,"sidata.contact_point is NaN!");
155 #endif
156         particles_in_contact = std::max(fi, fj) < 0.0;
157 
158         if(particles_in_contact) {
159           double contact_point_i[3], contact_point_j[3];
160           vectorSubtract3D(particle_j.gradient, particle_i.gradient, sidata.en);
161           vectorNormalize3D(sidata.en); //normalize
162 
163           double *const alpha_i = &sidata.contact_history[alpha1_offset];
164           double *const alpha_j = &sidata.contact_history[alpha2_offset];
165 
166           sidata.deltan = extended_overlap_algorithm(&particle_i, &particle_j, sidata.en, alpha_i, alpha_j,
167                 sidata.contact_point, contact_point_i, contact_point_j, sidata.delta);
168 
169           sidata.reff = sidata.radi*sidata.radj / (sidata.radi + sidata.radj);
170           if(curvatureLimitFactor > 0.0) {
171               int curvature_radius_method = meanCurvature ? 0 : 1;
172               double koefi = particle_i.calc_curvature_coefficient(curvature_radius_method, contact_point_i); //mean curvature coefficient
173               double koefj = particle_j.calc_curvature_coefficient(curvature_radius_method, contact_point_j); //mean curvature coefficient
174 #ifdef LIGGGHTS_DEBUG
175               if(std::isnan(koefi))
176                 error->one(FLERR,"sidata.koefi is NaN!");
177               if(std::isnan(koefj))
178                 error->one(FLERR,"sidata.koefj is NaN!");
179               if(std::isnan(sidata.reff))
180                 error->one(FLERR,"sidata.koefj is NaN!");
181 #endif
182               sidata.reff = get_effective_radius(sidata, atom->blockiness[iPart], atom->blockiness[jPart], koefi, koefj, curvatureLimitFactor, this->error);
183           }
184 #ifdef LIGGGHTS_DEBUG
185           if(std::isnan(vectorMag3D(contact_point_i)))
186             error->one(FLERR,"sidata.contact_point_i is NaN!");
187           if(std::isnan(vectorMag3D(contact_point_j)))
188             error->one(FLERR,"sidata.contact_point_j is NaN!");
189 
190           if(std::isnan(*alpha_i))
191             error->one(FLERR,"alpha_i is NaN!");
192           if(std::isnan(*alpha_j))
193             error->one(FLERR,"alpha_j is NaN!");
194           if(std::isnan(sidata.deltan))
195             error->one(FLERR,"sidata.deltan is NaN!");
196 #endif
197 
198           *particles_were_in_contact = SURFACES_INTERSECT;
199         } else
200           *particles_were_in_contact = SURFACES_CLOSE;
201        } else
202          *particles_were_in_contact = SURFACES_FAR;
203       *inequality_start = static_cast<double>(int_inequality_start);
204       return particles_in_contact;
205     }
206 
207     inline void surfacesIntersect(SurfacesIntersectData & sidata, ForceData&, ForceData&)
208     {
209       if(sidata.is_wall) {
210         sidata.reff = sidata.radi;
211         if(curvatureLimitFactor > 0.0) {
212           const int iPart = sidata.i;
213           particle_i.set(atom->x[iPart], atom->quaternion[iPart], atom->shape[iPart], atom->blockiness[iPart]);
214           int curvature_radius_method = meanCurvature ? 0 : 1;
215           double koefi = particle_i.calc_curvature_coefficient(curvature_radius_method, sidata.contact_point); //mean curvature coefficient
216           sidata.reff = get_effective_radius_wall(sidata, atom->blockiness[iPart], koefi, curvatureLimitFactor, this->error);
217         }
218       }
219       MathExtraLiggghtsNonspherical::surfacesIntersectNonSpherical(sidata, atom->x);
220     }
221 
222     inline void endSurfacesIntersect(SurfacesIntersectData&,TriMesh *, double * const) {}
223     inline void surfacesClose(SurfacesCloseData&, ForceData&, ForceData&){}
224     void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){}
225     void endPass(SurfacesIntersectData&, ForceData&, ForceData&){}
226     inline void tally_pp(double,int,int,int) {}
227     inline void tally_pw(double,int,int,int) {}
228 
229   protected:
230      double curvatureLimitFactor;
231      bool meanCurvature;
232      bool gaussianCurvature;
233   };
234 }
235 }
236 
237 #else  // SUPERQUADRIC_ACTIVE_FLAG
238 
239 // add dummy class
240 SURFACE_MODEL_DUMMY(SURFACE_SUPERQUADRIC, "Surface model used without support for superquadric particles")
241 
242 #endif // SUPERQUADRIC_ACTIVE_FLAG
243 
244 #endif // SURFACE_MODEL_SUPERQUADRIC_H_
245 #endif // SURFACE_MODEL
246