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 Alexander Podlozhnyuk 39 Copyright 2015- DCS Computing GmbH, Linz 40 ------------------------------------------------------------------------- */ 41 42 #ifndef LMP_SUPERQUADRIC_H 43 #define LMP_SUPERQUADRIC_H 44 45 #include <cstddef> 46 #include "math_extra.h" 47 48 struct Superquadric 49 { //interface structure for easy coding 50 double *center; 51 double *quat; 52 double *shape; 53 double *blockiness; 54 double gradient[3]; 55 double hessian[9]; 56 double rotation_matrix[9]; 57 double koef; 58 double shape_inv[3]; 59 bool isEllipsoid; 60 bool isCylinder; 61 bool useIntBlockiness; 62 63 void local2global(const double *input_coord, double *result); 64 void global2local(const double *input_coord, double *result); 65 void rotate_local2global(const double *input_coord, double *result); 66 void rotate_global2local(const double *input_coord, double *result); 67 68 void apply_homothety(double scale_factor, double *input_point, double *output_point); //shrink or blow a particle 69 double shape_function_local(const double *point); //calculates shape function value at the point, defined in the local(particle based) reference frame 70 double shape_function_global(const double *point); //calculates shape function value at the point, defined in the global reference frame 71 void shape_function_gradient_local(const double *input_coord, double *result); //calculates gradient of the shape function at the point, defined in the local(particle based) reference frame 72 void shape_function_gradient_global(const double *input_coord, double *result); //calculates gradient of the shape function at the point, defined in the global reference frame 73 void shape_function_hessian_local(const double *input_coord, double *result); //calculates hessian of the shape function at the point, defined in the local(particle based) reference frame 74 void shape_function_hessian_global(const double *input_coord, double *result); //calculates hessian of the shape function at the point, defined in the global reference frame 75 void map_point(const double *input_coord, double *result); // 76 void map_point_local(double *input_coord, double scale); 77 void reference_point(int iphi, int nphi, int itheta, int ntheta, double *point); 78 void pre_initial_estimate(const double *input_point, int nphi, int ntheta, int *iphi, int *itheta); 79 bool plane_intersection(const double *normal, const double *x0, double *result_point, double *point_of_lowest_potential); //finds the point of maximal penetration on the particle surface and the point of minimal shape function value on the surface 80 double line_intersection(const double *pointA, const double *pointB, double *closestPoint); //finds a point on edge incl.corners with a minimal value of particle shape function (point of lowest potential) 81 bool edge_intersection(const double *pointA, const double *pointB); //checks superquadric-edge intersection (only yes or no) 82 void tensor_quat_rotate(double *tensor, double *result); //rotate tensor 83 double calc_curvature_coefficient(int curvature_radius, const double *input_point); //calculates the mean curvature radius of a particle surface at a given point 84 void shape_function_props_local(const double *input_coord, double *f, double *grad, double *hess); //calculates shape function value, gradient and hessian at the point defined in local basis by one function 85 void shape_function_props_global(const double *input_coord, double *f, double *grad, double *hess); //calculates shape function value, gradient and hessian at the point defined in global basis by one function 86 87 void set_shape(double a, double b, double c); //sets particle shape parameters 88 void set_blockiness(double n1, double n2); //sets particle blockiness parameters 89 void calc_koef(); 90 91 double surface_line_intersection(bool use_alhpa, const double *start_point, const double *normal_vector, double alpha1, double *result); //calculates the intersection point between a line and particle surface with the Newton's method 92 double surface_line_intersection(const int max_num_iters, bool use_alhpa, const double *start_point, const double *direction_vector, double alpha1, double *result); 93 double surface_line_intersection1(const double *start_point, const double *direction_vector, double *result); 94 void ellipsoid_line_intersection(const double *start_point, const double *direction, double& alpha1, double& alpha2); 95 void ellipsoid_line_intersection_local(const double *start_point_local, const double *direction_vector_local, double& alpha1, double& alpha2); 96 double point_surface_projection(const int max_num_iters, bool use_alpha_ellipsoid, const double *start_point, double alpha1, double *result); 97 double calc_F(double *p, double *p0, double l, double *F); 98 SuperquadricSuperquadric99 Superquadric(): 100 center(NULL), 101 quat(NULL), 102 shape(NULL), 103 blockiness(NULL) 104 { 105 isEllipsoid = false; 106 isCylinder = false; 107 useIntBlockiness = false; 108 koef = 1.0; 109 } 110 Superquadric(double *center_, double *quat_, double *shape_, double *blockiness_); 111 void set(double *center_, double *quat_, double *shape_, double *blockiness_); 112 }; 113 114 #endif 115 116