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