1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
|
/*************************************************************************
* 2021 Jérôme Duriez, jerome.duriez@inrae.fr *
* 2023 Danny van der Haven, dannyvdhaven@gmail.com *
* This program is free software, see file LICENSE for details. *
*************************************************************************/
#ifdef YADE_LS_DEM
#pragma once
#include <lib/computational-geometry/MarchingCube.hpp>
#include <core/Shape.hpp>
#include <pkg/levelSet/RegularGrid.hpp>
namespace yade {
class LevelSet : public Shape {
private:
Vector3r center;
Real volume, lengthChar;
Vector3r inertia; // the eigenvalues of the inertia matrix: its diagonal expression in localAxes basis (in a Vector3r form here)
bool initDone;
bool initDoneMarchingCubes;
int nVoxInside;
vector<Vector3r> nodesOnRay; // computed in rayTrace
void assignSurfNodes(vector<Vector3r>);
vector<vector<vector<int>>> binarize(bool);
void postProcessNodes(); // applies necessary operations (e.g., computing sphericity) after surfNodes is filled
void init(); // compute nVoxInside, center, volume, and inertia
void rayTraceSurfNodes(const int&, const int&, const Real&); // fills surfNodes by ray tracing
void init_neighborsNodes(); // fills neighborsNodes
Real distanceInterpolation(const Vector3r&, const int&, const int&, const int&) const; // trilinear interpolation of distance in a given cell
bool rayTraceInCell(
const Vector3r&, const Vector3r&, const Vector3r&, const Vector3i&, const Real&); // handles the ray tracing from a given point in a given cell
vector<Vector3r> rayTrace(
const Vector3r&,
const Real&); // recursively calls rayTraceInCell, walking accross the whole grid along a ray starting from center and returning nodesOnRay
Real smearedHeaviside(Real) const;
Real volumeInsideThreshold(Real) const;
struct mcData { // Structure for holding marching cubes triangulation of level set particle
vector<Vector3r> triangles;
vector<Vector3r> normals;
int nbTriangles;
};
mcData marchingCubesData; // Actual marching cubes data holder
public:
Real minRad, maxRad; // for sphericity. Public for use in Ig2 for distant interactions
Real distance(const Vector3r&, const bool& unbound = false) const; // gives the distance from a point to the surface
Vector3r normal(const Vector3r&, const bool& unbound = false) const; // gives the outwards normal at some point
Real getVolume() override; // these 4 get*() may call init() if not already done, they can not be const-declared
Real getVolumeDep();
Vector3r getCenter();
Vector3r getInertia();
Real getSurface(Real epsilon = 1) const; // this one can be const-declared
Real getSurface_epsilon(Real) const; // handy function for getSurface
void computeMarchingCubes(); // Compute the marching cube triangulation for the LS shape
vector<Vector3r> getMarchingCubeTriangles(); // Retrieve marching cube triangles
vector<Vector3r> getMarchingCubeNormals(); // Retrieve marching cube normals
int getMarchingCubeNbTriangles(); // Retrieve marching cube number of triangles
virtual ~LevelSet() {};
// clang-format off
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(LevelSet,Shape,"A level set description of particle shape based on a :yref:`discrete distance field<LevelSet.distField>` (denoting $\\phi$ the underlying signed distance function) and :yref:`surface nodes<LevelSet.surfNodes>` [Duriez2021a]_ [Duriez2021b]_. See :ysrc:`examples/levelSet` for example scripts.",
((vector< vector< vector<Real> > >,distField,,Attr::readonly,"The signed (< 0 when inside) distance-to-surface function as a discrete scalar field on :yref:`lsGrid<LevelSet.lsGrid>`, with `distField[i][j][k]` corresponding to `lsGrid.gridPoint(i,j,k)`. From Python, slice this multi-dimensional list with care: while `distField[i][:][:]` corresponds to values on a x-cst plane, `distField[:][:][k]` is not at z-constant (use `[[distField[i][j][k] for j in ..] for i in ..]` instead)"))
((vector<Vector3r>,corners,,Attr::readonly,"The 8 corners of an axis-aligned bounding box, in local axes. It is computed once for all by :yref:`Bo1_LevelSet_Aabb` and used by the same Functor to get :yref:`Body.bound`."))
((vector<Vector3r>,surfNodes,,Attr::readonly,"Surface discretization in terms of (a list of) nodes, for the master-slave refined contact treatment in :yref:`Ig2_LevelSet_LevelSet_ScGeom`, previously coined boundNodes in [Duriez2021b]_. Expressed in local frame.")) // NB: just nodes as a name would "conflict" with many PFacet variables
((vector<vector<int>>,neighborsNodes,,Attr::readonly,"Denoting :yref:`n_neighborsNodes<LevelSet.n_neighborsNodes>` as $N$, neighborsNodes[i] is a $N+1$ elements iterable with i as the first item, then the indices of the $N$ closest neighbor nodes for every surface node i, ordered by increasing distance wrt to i."))
((int,n_neighborsNodes,-1,,"How many closest nodes we want to store for every one in :yref:`neighborsNodes<LevelSet.neighborsNodes>` (i.e., the $N$ therein). A negative value leaves the latter uncomputed i.e. empty"))
((Real,sphericity,-1,Attr::readonly,"Shape sphericity computed from boundary nodes and assuming both largest inscribed sphere and smallest circumscribed sphere have the origin (of local axes) as center."))
((bool,starLike,false,Attr::readonly,"Indicates whether :yref:`surface nodes<LevelSet.surfNodes>` are ray-traced considering that a ray can intersect the surface no more than once."))
((shared_ptr<RegularGrid>,lsGrid,new RegularGrid,Attr::readonly,"The :yref:`regular grid<RegularGrid>` carrying :yref:`distField<LevelSet.distField>`, in local axes."))
((bool,twoD,false,Attr::readonly,"True for z-invariant shapes. Serves to restrict the ray tracing of :yref:`surfNodes<LevelSet.surfNodes>` in the (x,y) plane."))
((Real,smearCoeff,1.5,,"Rules the smearing coefficient $\\varepsilon > 0$ of the Heaviside step function for a smooth integration of the particle's volume close to its surface (the higher $\\varepsilon$ the smoother, i.e. the more diffuse the surface in terms of volume integration). Given in reciprocal multiples of $R_{cell} = \\sqrt{3}/2 g_{grid}$ the half diagonal of the cells of the :yref:`lsGrid<LevelSet.lsGrid>` with $g_{grid}$ the :yref:`cell length<RegularGrid.spacing>`: $\\varepsilon = R_{cell}\\times 1/$ *smearCoeff* (smearing is deactivated if negative)."))
((bool,hasAABE,false,,"Flag to indicate whether an axis-aligned bounding ellipsoid (AABE) has been provided by the user. If true, you must specify :yref:`axisAABE<LevelSet.axisAABE>`. Only works for VLS-DEM."))
((Vector3r,axesAABE,Vector3r::Zero(),,"The half lengths of the principal axes of the axis-aligned bounding ellipsoid (AABE) of the level-set shape. Format (rx,ry,rz). Only works for VLS-DEM."))
,
minRad = std::numeric_limits<Real>::infinity();
maxRad = 0;
inertia = Vector3r::Zero();
initDone = false; // after hesitation, it is finally chosen to save the least of data, but to recall init() after a save/load
initDoneMarchingCubes = false;
lengthChar = -1;
volume = -1;
nVoxInside = -1;
center = Vector3r(std::numeric_limits<Real>::infinity(),std::numeric_limits<Real>::infinity(),std::numeric_limits<Real>::infinity());
createIndex(); // necessary for such a Shape-derived class, see https://yade-dem.org/doc/prog.html#indexing-dispatch-types
,
.def("binarize",&LevelSet::binarize,(boost::python::arg("surface")=false),R"""(Returns (for any user-defined purpose) a voxellised description of the shape volume. Selection (as being part of the particle) of voxels associated with a strictly zero value of distance depends on *surface*.
:param bool surface: whether to consider particle surface voxels (with a zero :yref:`distance<LevelSet.distField>` value) as being part of the volume (if True) or not (if False)
:returns: a 3D grid of 0/1 integers whose [i][j][k] term is 1, resp. 0, when the voxel associated to :yref:`lsGrid.gridPoint(i,j,k)<RegularGrid.gridPoint>` (see :yref:`volume<LevelSet.volume>` for an exact definition of that association) is inside, resp. outside, the particle. Modulo *surface* considerations, the returned 0/1 values actually correspond to $(\phi-|\phi|)/(2 \phi)$
)""")
.def("volume",&LevelSet::getVolumeDep,"A 2025/09-deprecated version of :yref:`getVolume<LevelSet.getVolume>`. Please use the latter instead.")
.def("getVolume",&LevelSet::getVolume,"Returns the shape volume as the negative domain of the :yref:`level set function<LevelSet.distField>`, in a voxellised fashion, where a negative value of the level set, :yref:`distField[i][j][k]<LevelSet.distField>` $\\leq$ 0, is considered to correspond to a material cubic voxel of side :yref:`lsGrid.spacing<RegularGrid.spacing>` and centered at :yref:`lsGrid.gridPoint(i,j,k)<RegularGrid.gridPoint>`. Smearing considerations may apply as per :yref:`smearCoeff<LevelSet.smearCoeff>`.")
.def("center",&LevelSet::getCenter,"The center of mass of the :yref:`volume<LevelSet.volume>` (considering obviously an uniform density for this volume), in local axes (for verification purposes, by comparison with the origin).")
.def("inertia",&LevelSet::getInertia,"The diagonal coefficients (i.e., eigenvalues, in a consistent workflow) of the geometric inertia matrix (the one considering the infinitesimal volume as the integrand, instead of infinitesimal mass) of the particle :yref:`volume<LevelSet.getVolume>`, as a (xx,yy,zz) Vector3r.")
// .def("nodesInCell",&LevelSet::getNodesInCellCube,(boost::python::args("i", "j", "k")),"Which boundary nodes belong to a given grid cube (given by its i,j,k indices)")
.def("distance",&LevelSet::distance,(boost::python::arg("pt"),boost::python::arg("unbound")=false),"Distance to surface at *pt*, with pt being expressed in the local frame. The 'unbound' flag (if True) allows the computation of distance values outside of the :yref:`grid<LevelSet.lsGrid>` extents, which otherwise returns NaN together with an error.")
.def("normal",&LevelSet::normal,(boost::python::arg("pt"),boost::python::arg("unbound")=false),"Unit normal vector to the surface, at some pt. Local frame applies to both output normal and input pt. Has an 'unbound' flag signaling whether to allow (if True) the computation of the normal outside of the :yref:`grid<LevelSet.lsGrid>` extents.")
.def("assignSurfNodes",&LevelSet::assignSurfNodes,(boost::python::arg("nodes")),"Assigns given *nodes* (as a list of Vector3r points conforming local frame) to :yref:`surfNodes<LevelSet.surfNodes>`, erasing pre-existing ones if any.")
.def("rayTrace",&LevelSet::rayTrace,(boost::python::arg("ray"),boost::python::arg("nodesTol")=50),"Performs one ray tracing along *ray* (to be given in cartesian coordinates, see spher2cart() if a conversion is useful from spherical ones; norm should not matter) and returns the obtained surface points (without modifying :yref:`surfNodes<LevelSet.surfNodes>`). Argument *nodesTol* has the same meaning as in :yref:`rayTraceSurfNodes<LevelSet.rayTraceSurfNodes>`, see therein. Provided for debugging or shape analysis purposes.")
.def("rayTraceSurfNodes",&LevelSet::rayTraceSurfNodes, (boost::python::arg("nSurfNodes")=102, boost::python::arg("nodesPath")=2, boost::python::arg("nodesTol") = 50),R"""(Defines all surface nodes (i.e., :yref:`surfNodes<LevelSet.surfNodes>`) by ray tracing, erasing pre-existing ones if any, according to given arguments
:param int nSurfNodes: how many rays are to be launched. Please use a perfect square + 2 if not :yref:`twoD<LevelSet.twoD>` and if *nodesPath* = 1. Previously coined nNodes in [Duriez2021b]_
:param int nodesPath: defines how the space of spherical coordinates $(\theta \in [0;\pi] ,\varphi\in [0;2 \pi])$ is discretized when ray tracing the boundary nodes: 1 gives a rectangular partition of that space, plus two nodes at $\theta = 0 [\pi]$; 2 locates the nodes along a spiral path [Duriez2021a]_
:param real nodesTol: tolerance coefficient for accepting (if $|\phi| / L <$ nodesTol $\times$ numeric precision with $\phi$ the return value of :yref:`distance<LevelSet.distance>` and $L$ a body-characteristic length taken as $\sqrt[3]{V}$ with $V$ the :yref:`volume<LevelSet.volume>`, or $\sqrt{V/g_{grid}}$ with $g_{grid}$ the grid :yref:`spacing<RegularGrid.spacing>` if :yref:`twoD<LevelSet.twoD>`) surface nodes proposed by the ray tracing algorithm
)""")
.def("getSurface",&LevelSet::getSurface,(boost::python::arg("epsilon")=1.),"Compute particle surface by differentiation of Level Set volume around $\\phi = 0$, where best results have usually been obtained using :yref:`smearCoeff<LevelSet.smearCoeff>` = 1. Given argument *epsilon* (as a length quantity) serves to initialize the recursive limit search in the derivative expression (see [Duriez2025]_)")
.def("computeMarchingCubes",&LevelSet::computeMarchingCubes,"Compute or recompute the triangulation of the particle surface after using the Marching Cubes algorithm on :yref:`distField<LevelSet.distField>`.")
.def("marchingCubesVertices",&LevelSet::getMarchingCubeTriangles,"Returns the vertices for a surface triangulation obtained after executing the Marching Cubes algorithm on :yref:`distField<LevelSet.distField>`.")
.def("marchingCubesNormals",&LevelSet::getMarchingCubeNormals,"Returns the normals for a surface triangulation obtained after executing the Marching Cubes algorithm on :yref:`distField<LevelSet.distField>`.")
.def("marchingCubesNbTriangles",&LevelSet::getMarchingCubeNbTriangles,"Returns the number of triangles forming the surface triangulation as per the Marching Cubes algorithm (executed on :yref:`distField<LevelSet.distField>`).")
)
// clang-format on
REGISTER_CLASS_INDEX(LevelSet, Shape); // necessary for such a Shape-derived class, see https://yade-dem.org/doc/prog.html#indexing-dispatch-types
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(LevelSet);
} // namespace yade
#endif // YADE_LS_DEM
|
