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