[mpm] Add ParticleSorter

[xuchenhan-tri]

In addition:

SpGrid reports the flags it uses.
Add the missing set_zero() method for GridData.

---

This change is 

[xuchenhan-tri]

+(release notes: none)
+@amcastro-tri for feature review please. This PR is getting close to touch the core transfer algorithm, and I figure you might be interested in understanding how it works.

Reviewable status: LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

[amcastro-tri]

Working... I left a couple of questions.
Mostly I'd need some variables definitions in order to dig into the implementation of Sort().

Reviewed 5 of 8 files at r1, all commit messages.
Reviewable status: 8 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

---

multibody/mpm/spgrid.h line 28 at r1 (raw file):

/* A subset of SPGrid flags and configs we care about. */
struct SpGridFlags {

per this PR, this is struct and accessor below are not used.
Remove.

---

multibody/mpm/particle_sorter.h line 25 at r1 (raw file):

};

using Ranges = std::vector<Range>;

nit, I'd personally advocate for a more explicit name like say RangesVector to make sure it differentiates by more than one character from Range and also it comunicates quickly there's a std::vector underneath.
I'd also be okay with using std::vector<Range> all around.

---

multibody/mpm/particle_sorter.h line 32 at r1 (raw file):

 @pre data is sorted in ascending order and has at least two elements.
 @pre ranges != nullptr */
void ConvertToRanges(const std::vector<int>& data, Ranges* ranges);

can the input vector have repeated elements? say {1, 3, 3, 6, 9}? would that lead to an "empty" range? is that allowed?
Simply document requirements and/or add checks if possible.

---

multibody/mpm/particle_sorter.h line 38 at r1 (raw file):

 Given a set of MPM particles and the background grid, recall that we define the
 base node of a particle to be the grid node that is closest to the particle

Would this be correct?

Suggestion:

 Given a set of MPM particles and the background grid, recall that we define the
 base node of a particle to be the grid node that is closest (in the max norm) to the particle

---

multibody/mpm/particle_sorter.h line 59 at r1 (raw file):

  template <typename T, typename SpGrid>
  void Sort(const SpGrid& spgrid, double dx,
            const std::vector<Vector3<T>>& particle_positions) {

not sure if documented elsewhere, but could you say or reference what's the phisical space occupied by an SpGrid with grid size dx?
Do the grid's coordinates always go from 0 to kMaxGridSize*dx? i.e. the origin is always at the corner of the grid?

In other words, I guess the SpGrid defines a frame (definition?) and the particle_poisions are expressed in that frame?

---

multibody/mpm/particle_sorter.h line 216 at r1 (raw file):

  std::vector<int> sentinel_particles_;
  std::vector<int> data_indices_;
  std::vector<uint64_t> base_node_offsets_;

maybe I could guess from the impl, but I think a good idea to document here locally the intent of these individual variables. Is the comment at the top supposed to apply to sentinel_particles_ or to all of these three arrays?

Also,w hat's the difference between an "index" (as stored in data_indices_) and an "offset" (as stored by base_node_offsets_).
I know what an offset is from the docs of SpGrid, but I don't know what an "index" is.

Could you also provide a local definitiion for "sentinel particles"? or maybe reference the paper that provides that definition?

It'd help a lot with the review and I believe also with the future maintenance of this code.

---

multibody/mpm/particle_sorter.h line 220 at r1 (raw file):

   index of the particle. */
  std::vector<uint64_t> particle_sorters_;
  std::array<Ranges, 8> colored_ranges_;

ditto for these, probably local comments that explain their specific purpose would serve best those of us unfamiliar with this code, rather than a generic header for a group of variables.

---

multibody/mpm/grid_data.h line 116 at r1 (raw file):

  /* Sets all floating point data to zero, and set the index to be inactive.*/
  void set_zero() {

without reading the docs this methods sounds a lot like a harmless Eigen setZero().
IMO, the most "surprising" or really releveant effect of this function is to set to inactive, with the side effect of setting things to zero.

In other words, if the method is called set_inactive() I will probably be more prone to reading the docs rather than guessing this function's behavior (erroneously as just setting things to zero without changing index state).

Another option is to be very verbose and call this set_inactive_and_zero(), but I believe set_inactive() should suffice.

[amcastro-tri]

Reviewable status: 7 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

---

multibody/mpm/spgrid.h line 28 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

per this PR, this is struct and accessor below are not used.
Remove.

nvm

[xuchenhan-tri]

Reviewable status: 6 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

---

multibody/mpm/grid_data.h line 116 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

without reading the docs this methods sounds a lot like a harmless Eigen setZero().
IMO, the most "surprising" or really releveant effect of this function is to set to inactive, with the side effect of setting things to zero.

In other words, if the method is called set_inactive() I will probably be more prone to reading the docs rather than guessing this function's behavior (erroneously as just setting things to zero without changing index state).

Another option is to be very verbose and call this set_inactive_and_zero(), but I believe set_inactive() should suffice.

There are two reasons that this is named set_zero().

1. It's required by (and called from) the SpGrid class. See class documentation of SpGrid in spgrid.h.
2. It's important all floating point values are zero, because new values will later be added to this struct without resetting things to zero.

If you think the naming is confusing, I can be persuaded to rename this to reset and remove the existing reset function.

---

multibody/mpm/particle_sorter.h line 25 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

nit, I'd personally advocate for a more explicit name like say RangesVector to make sure it differentiates by more than one character from Range and also it comunicates quickly there's a std::vector underneath.
I'd also be okay with using std::vector<Range> all around.

Done

---

multibody/mpm/particle_sorter.h line 32 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

can the input vector have repeated elements? say {1, 3, 3, 6, 9}? would that lead to an "empty" range? is that allowed?
Simply document requirements and/or add checks if possible.

Done

---

multibody/mpm/particle_sorter.h line 38 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

Would this be correct?

That would be correct, but I don't see the need to specify the norm. The result is the same with any reasonable norm (e.g. euclidean norm).

---

multibody/mpm/particle_sorter.h line 59 at r1 (raw file):
SpGrid doesn't have the notion of dx. It's a lattice of integer valued 3 vectors.

Do the grid's coordinates always go from 0 to kMaxGridSize*dx? i.e. the origin is always at the corner of the grid?

This is abstracted away (as much as we could). As the user of SpGrid, you can basically treat it as an infinite lattice of Vector3i with a bijection from Vector3i to uint64t.

In other words, I guess the SpGrid defines a frame (definition?) and the particle_poisions are expressed in that frame?

I clarified the frame for particle positions.

---

multibody/mpm/particle_sorter.h line 216 at r1 (raw file):
All private member variables without documentation is explained in the accessor. I've made that clear.

Is the comment at the top supposed to apply to sentinel_particles_ or to all of these three arrays?

Just sentinel_particles.

Could you also provide a local definitiion for "sentinel particles"? or maybe reference the paper that provides that definition?

No paper discusses this concept AFAIK. This is an implementation detail. I made an attempt to clarify what a sentinel particle a bit more.

---

multibody/mpm/particle_sorter.h line 220 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

ditto for these, probably local comments that explain their specific purpose would serve best those of us unfamiliar with this code, rather than a generic header for a group of variables.

This is already explained in the accessor().

[amcastro-tri]

Reviewable status: 7 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

---

multibody/mpm/grid_data.h line 116 at r1 (raw file):

1. It's required by (and called from) the SpGrid class.

I know, but that's under our control, you can always change the function's name both here and in spgrid.h.

2. It's important all floating point values are zero

Yeah, I do see that also.

If you think the naming is confusing, I can be persuaded to rename this to reset

Well, you have a reset() now that sets things to NaN, I imagine for good reasons.

I think you could keep your reset() and have this method to have a more explicit funtion name. Something like zero_and_deactivate() could do?
I just think that within local scope (say in spgrid.h), saying set_zero() gloses over the important detail that now this data will be marked inactive.

You could also consider two separate calls if you hate the compound name, i.e. a set_zero() (sets all floating points to zero, leaves index untouched) and a set_inactive() (marks the index/data inactive).

---

multibody/mpm/particle_sorter.h line 38 at r1 (raw file):

Previously, xuchenhan-tri wrote…

That would be correct, but I don't see the need to specify the norm. The result is the same with any reasonable norm (e.g. euclidean norm).

ah, you are right. This is good. I extrapolated my mental picture of a box around the base node to a box associated with the max norm.

---

multibody/mpm/particle_sorter.h line 59 at r1 (raw file):

As the user of SpGrid, you can basically treat it as an infinite lattice of Vector3i

Ah, are you sayint that these i-j-k indexes can be negative also? and thus SpGrid mathematically spans all of R^3?

---

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):

The order in which the particle data is preferrablly accessed.

Tryint to reconcilliate this with the doc for base_node_offsets().

Would it be true to say:

Returns a vector toaid the recursion of particle data in "sorted" order. This is the in which data is preferrably accessed when working with SpGrid.

I guess what I liked to get from this comment (mine might be wrong) is:

1. "Order", is this the "sorted order" referenced in base_node_offsets()?
2. "Prefered", is this a requirement for the best usage of SpGrid?

---

multibody/mpm/particle_sorter.h line 66 at r2 (raw file):

  template <typename T, typename SpGrid>
  void Sort(const SpGrid& spgrid, double dx,
            const std::vector<Vector3<T>>& particle_positions) {

Tryint to reconcilliate how indexes returned by data_indices() are assigned.
Are those indexes assigned based on the order in which particles are specified by particle_positions? i.e. the first position has index 0, the second index 1, etc.?

---

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):

  std::vector<uint64_t> GetActiveBlockOffsets() const;

  /* The order in which the particle data is preferrablly accessed.

minor,

Suggestion:

preferablly

---

multibody/mpm/particle_sorter.h line 194 at r2 (raw file):

    for (int p = 0; p < num_particles; ++p) {
      const auto& particle_data = foo_data[data_indices[p]];
      // do something with particle_data.

ah, this is great. I guess what I didn't get from SpGrid's docs is that the offsets are not contiguous then? I was thinking I'd use offsets as 1D indexes. If true, maybe worth pointing out either here on SpGrid's docs.

---

multibody/mpm/particle_sorter.h line 209 at r2 (raw file):

   kernels, without write hazards between concurrent threads.

   Specifically, for each i from 0 to 7, colored_ranges()[i] is a collection of

I love this second paragraph, very clear.

Could you say something about this fixed (8) number of colors? why couldn't this be a different number? do we always have 8 colors and possibly empty ranges? or could we ever have more than 8 colors?
If needed, cite reference if you have it.

[xuchenhan-tri]

Reviewable status: 7 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

---

multibody/mpm/grid_data.h line 116 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

1. It's required by (and called from) the SpGrid class.

I know, but that's under our control, you can always change the function's name both here and in spgrid.h.

2. It's important all floating point values are zero

Yeah, I do see that also.

If you think the naming is confusing, I can be persuaded to rename this to reset

Well, you have a reset() now that sets things to NaN, I imagine for good reasons.

I think you could keep your reset() and have this method to have a more explicit funtion name. Something like zero_and_deactivate() could do?
I just think that within local scope (say in spgrid.h), saying set_zero() gloses over the important detail that now this data will be marked inactive.

You could also consider two separate calls if you hate the compound name, i.e. a set_zero() (sets all floating points to zero, leaves index untouched) and a set_inactive() (marks the index/data inactive).

I could live without the current reset function and have the default state of the grid state to be zero and inactive. In that case reset is properly named.

---

multibody/mpm/particle_sorter.h line 59 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

As the user of SpGrid, you can basically treat it as an infinite lattice of Vector3i

Ah, are you sayint that these i-j-k indexes can be negative also? and thus SpGrid mathematically spans all of R^3?

Not technically R^3 because they are integer valued, but yes, they can be negative

---

multibody/mpm/particle_sorter.h line 66 at r2 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

Tryint to reconcilliate how indexes returned by data_indices() are assigned.
Are those indexes assigned based on the order in which particles are specified by particle_positions? i.e. the first position has index 0, the second index 1, etc.?

I believe this resolved per our f2f, but feel free to continue the discussion here if you disagree.

---

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):
Working.

• "Order", is this the "sorted order" referenced in base_node_offsets()?

Yes

• "Prefered", is this a requirement for the best usage of SpGrid?

It's not neither a requirement nor necessarily absolutely the the "best", but it's the best I (and many before myself) can come up with.

---

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

minor,

Working

---

multibody/mpm/particle_sorter.h line 194 at r2 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

ah, this is great. I guess what I didn't get from SpGrid's docs is that the offsets are not contiguous then? I was thinking I'd use offsets as 1D indexes. If true, maybe worth pointing out either here on SpGrid's docs.

I think this is resolved from our f2f? Feel free to reopen and request changes if you disagree.

---

multibody/mpm/particle_sorter.h line 209 at r2 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

I love this second paragraph, very clear.

Could you say something about this fixed (8) number of colors? why couldn't this be a different number? do we always have 8 colors and possibly empty ranges? or could we ever have more than 8 colors?
If needed, cite reference if you have it.

Working