geometry: Use field gradient and surface normal to filter ContactSurface. #12730
Conversation
DamrongGuoy
force-pushed
the
field_gradient
branch
2 times, most recently
from
94b0efd
to
c8cc042
Compare
There was a problem hiding this comment.
+@SeanCurtis-TRI for feature review please.
CC: @sherm1 , @tehbelinda
If you create geometry's doxygen bazel-bin/doc/doxygen --quick geometry, you should see this Module doc:
- Drake
- Modules
- Geometric Representations
- Proximity Queries
- Collision Queries
- Contact Surface
It explains how thin objects can cause problems.
- Contact Surface
- Collision Queries
- Proximity Queries
- Geometric Representations
- Modules
Reviewable status: LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @SeanCurtis-TRI)
There was a problem hiding this comment.
I'm running off to CPR training...here are some comments to resolve.
Reviewed 2 of 24 files at r1.
Reviewable status: 5 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy and @SeanCurtis-TRI)
geometry/proximity/mesh_field.h, line 58 at r2 (raw file):
@param e The index of the element. @param b The barycentric coordinates. @note Currently the abstract class MeshField has only one concrete subclass
I don't think this should go in MeshField.
- Since we originally authored this, it has become blazingly clear that we have no expectations for anything but linear approximations of the pressure field.
- The only concrete implementation of this interface ignores part of the interface that suggests as a "generalization", it's not very applicable.
- I'm almost to the point where I think we should simply kill
MeshFieldand just run withMeshFieldLinearalone. However, definitely not in this PR.
So, for all those reasons, I believe we should simply leave this API alone in this PR.
geometry/proximity/mesh_field_linear.h, line 23 at r2 (raw file):
At the high level, I'd say something like this (before you get into all the details below):
In other words, it represents a smooth, continuous field as a continuous piecewise linear function with a discontinuous gradient; the field value changes linearly within each element and the gradient is constant within each element.
This calls out the fact that the function is continuous and strongly implies it isn't smooth. I want all of the interesting, relevant properties captured in the first summary paragraph and all of the details can simply elaborate on those properties.
geometry/proximity/mesh_field_linear.h, line 32 at r2 (raw file):
For a linear triangle or tetrahedron in 3-D, we use barycentric coordinates bᵢ, where bᵢ is the weight of the iᵗʰ vertex, and their summation is 1. The barycentric coordinates B is:
nit s/coordinates/coordinate (singular)
geometry/proximity/mesh_field_linear.h, line 41 at r2 (raw file):
I'd like to discuss the notation you have here -- we can rope @mitiguy in on this one. But I think you have a problem with the notation that prevents clarity. Here's the gist of my issue.
- We have a point Q.
- It lies "in" a simplicial element, E.
- We can describe it as a cartesian vector: p_FQ_F (w.r.t. some frame F)
- We can also describe it as a barycentric vector: b_EQ (w.r.t. element E)
The key to this is a point is a point. And a barycentric coordinate is analogous to a cartesian coordinate. Neither of them are the point itself, but merely a way of describing it which depends on conventions and reference frames. So, with that in mind, the documentation above would become:
For a point Q which lies "in" the simpicial element E with barycentric
coordinates b_EQ, Q can be expressed in cartesian coordinates as:
p_MQ := ∑b_EQᵢ*r_MVᵢwhere `r_MVᵢ is the position vector of vertex Vᵢ expressed in Frame M.
This cleans up the fact that r_M isn't part of our monogram notation and B isn't a point, but a coordinate.
geometry/proximity/mesh_field_linear.h, line 47 at r2 (raw file):
<h3> Field Approximation </h3> At the point r_M corresponding to the barycentric coordinates B, the weighted
See point above about r_M and B.
There was a problem hiding this comment.
Reviewable status: 5 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy and @SeanCurtis-TRI)
geometry/proximity/mesh_field_linear.h, line 41 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
I'd like to discuss the notation you have here -- we can rope @mitiguy in on this one. But I think you have a problem with the notation that prevents clarity. Here's the gist of my issue.
- We have a point Q.
- It lies "in" a simplicial element, E.
- We can describe it as a cartesian vector: p_FQ_F (w.r.t. some frame F)
- We can also describe it as a barycentric vector: b_EQ (w.r.t. element E)
The key to this is a point is a point. And a barycentric coordinate is analogous to a cartesian coordinate. Neither of them are the point itself, but merely a way of describing it which depends on conventions and reference frames. So, with that in mind, the documentation above would become:
For a point Q which lies "in" the simpicial element E with barycentric
coordinates b_EQ, Q can be expressed in cartesian coordinates as:
p_MQ := ∑b_EQᵢ*r_MVᵢwhere `r_MVᵢ is the position vector of vertex Vᵢ expressed in Frame M.
This cleans up the fact that
r_Misn't part of our monogram notation andBisn't a point, but a coordinate.
BTW I'd happily swap b_EQ for just b_Q (as I believe that's notation that I've used elsewhere (e.g., surface_mesh.h, and most other places).
There was a problem hiding this comment.
Another breakpoint.
Reviewed 1 of 24 files at r1.
Reviewable status: 12 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy and @SeanCurtis-TRI)
geometry/proximity/mesh_field_linear.h, line 21 at r2 (raw file):
/** %MeshFieldLinear represents a field variable defined on a simplicial
BTW We have a very concrete notational challenge to work out. And that is the index of a quantity in the set and in the face. So, a mesh has the ordered set of vertices V = {v0, v1, v2..., vN-1}. Thus we could refer to vertex v2. However, an element has "vertices" (indices into the set of vertices) (i0, i1, ...). THen, finally, we have the vertices the faces refers to as something like V(i0), V(i1), V(i2). Such that V(i0) is not necessarily v0. This is not necessarily the notation I'd advocate but I think we need to be clear and consistent about distinguishing those.
I could also accept the idea of saying the actual vertices in the set are V0, V1, V2, ..... and the v0 = V(i0) such that v0 is an alias for some vertex Vk. I'm open to any number of things. But I'm not particularly pleased with the current blurred line.
--- As per our f2f discussion -- we've reached consensus on the problem statement and recognize that resolution is no trivial. So, I've copied my notes below and will rely on your judgement as to an appropriate solution.
f(x) x \in domain of mesh
f*(x) -- linear approximation of f
- Mesh
- V = {v0, v1, ...} - ordered set of vertices
- E = {{i0, i1, i2}, ... } - ordered set of triples -- indices into vertices
- b(x) --> barycentric coordinates of the point X (assuming X in mesh's domain)
- element(x) --> defines element
- f_samples = U = {f(vi) for vi \in V}
- grad_f_samples = gradU = { ?(e_i) for e in E}
such that
f*(x) = b0(x) * U(i0) + b1(x) * U(i1) + b2(x) * U(i2)
gradf*(x) = gradU(element(x))
geometry/proximity/mesh_field_linear.h, line 48 at r2 (raw file):
At the point r_M corresponding to the barycentric coordinates B, the weighted sum of the field values uᵢ's at vertex Vᵢ's is the piecewise linear
Probably worth defining u up when you talk about "store one field value per vertex".
geometry/proximity/mesh_field_linear.h, line 49 at r2 (raw file):
nit: this sentence is overly convoluted. And because of that, it lost a bit at the end:
...of the field u on the element E at the point P.
(With me opting to call the point P based on my earlier comments. It's neither point r, M, B, u, or E, which are the only symbols available in the current text.)
geometry/proximity/mesh_field_linear.h, line 53 at r2 (raw file):
uᵉ = b₀ * u₀ + b₁ * u₁ + b₂ * u₂ for triangle E, = b₀ * u₀ + b₁ * u₁ + b₂ * u₂ + b₃ * u₃ for tetrahedron E.
BTW I'd advocate not spelling out each variation for tri and tet. I think it's worth noting that the dimension of a barycentric coordinate is 3 and 4 (for tri and tet, respectively) and then express all the math simply as a sum over ∑b_Qᵢ.
geometry/proximity/mesh_field_linear.h, line 56 at r2 (raw file):
<h3> Gradient </h3> Consider each barycentric coordinate bᵢ as a linear scalar field on an
BTW The more I read about this documentation, the more I wonder what it's value is. There's a lot of detail that doesn't empower the user of this class. What they need to know:
- This is a piecewise linear function
- It has a piecwise constant gradient.
- The domain of each piece is an element.
- The function value at a point P in an element is queried via the barycentric coordinates of point b_P in that element. (And given r_MP, where do they get b_P)?
- There's a cheap backdoor for evaluating the function at vertex values.
That's essentially everything they need to know.
There's a bit more detail in the constructor (i.e., sample f(x) at each vertex and pass it in), but other than that. It seems this documentation should be quite simple.
Per our f2f -- let's
geometry/proximity/mesh_field_linear.h, line 64 at r2 (raw file):
Note that ∇bᵢ is constant on E and depends on the shape of the triangle or tetrahedron E.
BTW This note is factual and true, but not complete, which might lead to apparent untruth. It depends on the shape of the element and the values of u at the element's vertices.
geometry/proximity/mesh_field_linear.h, line 121 at r2 (raw file):
std::vector<FieldValue>& mutable_values() { return values_; } /**
nit: This shouldn't be public.
geometry/proximity/mesh_field_linear.h, line 136 at r2 (raw file):
@param field The field for comparison. @returns `true` if the given field is equal. @note Requires `MeshField field` to be MeshFieldLinear.
nit: Please put in a TODO on this changing the parameter type to be MeshFieldLinear. As long as it's a requirement, we should have it check at compile time.
Or is there a particular reason why it's expressed in terms of the base class?
geometry/proximity/mesh_field_linear.h, line 180 at r2 (raw file):
template <class FieldValue, class MeshType> void MeshFieldLinear<FieldValue, MeshType>::CalcGradientField() {
Any particular reason these are not defined inline when all of the other functions are?
geometry/proximity/mesh_intersection.h, line 479 at r2 (raw file):
thin rigid plate N ┌─┐
BTW I had a passing thought. If you want to emphasize the contact surface, you could do it something like this:
┌┄┐
┊ ┊ soft ball M
┊ ┊ ● ● ● ●
┊ ║● ●
⇦┃↘║⇨ ●
●⇦┃↘║↘ ●
● ⇦┃ ║⇨ ↘ ●
● ⇦┃ ║⇨ ↘ ●
● ⇦┃→║⇨ → → → ●
● ⇦┃ ║⇨ ↗ ● ↗ pressure gradient ∇p_M in M
● ⇦┃ ║⇨ ↗ ● ⇨ surface normal f_N on N
●⇦┃↗║↗ ● ║ suitable intersecting surface
⇦┃↗║⇨ ● ┃ unsuitable intersecting surface
┊ ║● ● ┊ non-intersecting surface
┊ ┊ ● ● ● ●
┊ ┊
└┄┘
geometry/proximity/mesh_intersection.h, line 496 at r2 (raw file):
│ │ └─┘
BTW These illustrations are gorgeous. Two thumbs way up!
geometry/proximity/mesh_intersection.h, line 504 at r2 (raw file):
However, there is no single angle threshold that works for all cases. For example, a rigid box N penetrates into a soft ball M in the following picture has triangles on its left side and right side that have
nit: grammar problem "a rigid box ... penetrates ... has triangles..." Two verbs on the box, but no conjunction. Just adding "and" would resolve that: "...following picture and has triangles..."
There was a problem hiding this comment.
Reviewable status: 7 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @DamrongGuoy and @SeanCurtis-TRI)
geometry/proximity/mesh_field_linear.h, line 21 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW We have a very concrete notational challenge to work out. And that is the index of a quantity in the set and in the face. So, a mesh has the ordered set of vertices V = {v0, v1, v2..., vN-1}. Thus we could refer to vertex v2. However, an element has "vertices" (indices into the set of vertices) (i0, i1, ...). THen, finally, we have the vertices the faces refers to as something like V(i0), V(i1), V(i2). Such that V(i0) is not necessarily v0. This is not necessarily the notation I'd advocate but I think we need to be clear and consistent about distinguishing those.
I could also accept the idea of saying the actual vertices in the set are V0, V1, V2, ..... and the v0 = V(i0) such that v0 is an alias for some vertex Vk. I'm open to any number of things. But I'm not particularly pleased with the current blurred line.
--- As per our f2f discussion -- we've reached consensus on the problem statement and recognize that resolution is no trivial. So, I've copied my notes below and will rely on your judgement as to an appropriate solution.
f(x) x \in domain of mesh f*(x) -- linear approximation of f - Mesh - V = {v0, v1, ...} - ordered set of vertices - E = {{i0, i1, i2}, ... } - ordered set of triples -- indices into vertices - b(x) --> barycentric coordinates of the point X (assuming X in mesh's domain) - element(x) --> defines element - f_samples = U = {f(vi) for vi \in V} - grad_f_samples = gradU = { ?(e_i) for e in E} such that f*(x) = b0(x) * U(i0) + b1(x) * U(i1) + b2(x) * U(i2) gradf*(x) = gradU(element(x))
Done. I hope the notation is better now.
geometry/proximity/mesh_field_linear.h, line 23 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
At the high level, I'd say something like this (before you get into all the details below):
In other words, it represents a smooth, continuous field as a continuous piecewise linear function with a discontinuous gradient; the field value changes linearly within each element and the gradient is constant within each element.
This calls out the fact that the function is continuous and strongly implies it isn't smooth. I want all of the interesting, relevant properties captured in the first summary paragraph and all of the details can simply elaborate on those properties.
Done.
geometry/proximity/mesh_field_linear.h, line 47 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
See point above about
r_MandB.
Done.
geometry/proximity/mesh_field_linear.h, line 48 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
Probably worth defining
uup when you talk about "store one field value per vertex".
Done. Also changed u to f per f2f discussion.
geometry/proximity/mesh_field_linear.h, line 53 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW I'd advocate not spelling out each variation for tri and tet. I think it's worth noting that the dimension of a barycentric coordinate is 3 and 4 (for tri and tet, respectively) and then express all the math simply as a sum over
∑b_Qᵢ.
Per f2f discussion, now I try b₀(Q)F[i₀] + b₁(Q)F[i₁] + b₂(Q)F[i₂]. The ∑ is hard due to "index of index: i₀, i₁, i₂" notation.
geometry/proximity/mesh_field_linear.h, line 64 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW This note is factual and true, but not complete, which might lead to apparent untruth. It depends on the shape of the element and the values of u at the element's vertices.
∇uᵉ depends on the values of u at vertices. ∇bᵢ does not depend on values of u.
geometry/proximity/mesh_field_linear.h, line 180 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
Any particular reason these are not defined inline when all of the other functions are?
Done. Changed to inline. I moved this function back and forth between public and private. Now it's public, so it's inline now.
geometry/proximity/mesh_intersection.h, line 479 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW I had a passing thought. If you want to emphasize the contact surface, you could do it something like this:
┌┄┐ ┊ ┊ soft ball M ┊ ┊ ● ● ● ● ┊ ║● ● ⇦┃↘║⇨ ● ●⇦┃↘║↘ ● ● ⇦┃ ║⇨ ↘ ● ● ⇦┃ ║⇨ ↘ ● ● ⇦┃→║⇨ → → → ● ● ⇦┃ ║⇨ ↗ ● ↗ pressure gradient ∇p_M in M ● ⇦┃ ║⇨ ↗ ● ⇨ surface normal f_N on N ●⇦┃↗║↗ ● ║ suitable intersecting surface ⇦┃↗║⇨ ● ┃ unsuitable intersecting surface ┊ ║● ● ┊ non-intersecting surface ┊ ┊ ● ● ● ● ┊ ┊ └┄┘
Done. Thanks.
DamrongGuoy
force-pushed
the
field_gradient
branch
from
bd1d7d0
to
59677b9
Compare
There was a problem hiding this comment.
Reviewable status: 7 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy, @mitiguy, and @SeanCurtis-TRI)
geometry/proximity/mesh_field.h, line 58 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
I don't think this should go in
MeshField.
- Since we originally authored this, it has become blazingly clear that we have no expectations for anything but linear approximations of the pressure field.
- The only concrete implementation of this interface ignores part of the interface that suggests as a "generalization", it's not very applicable.
- I'm almost to the point where I think we should simply kill
MeshFieldand just run withMeshFieldLinearalone. However, definitely not in this PR.So, for all those reasons, I believe we should simply leave this API alone in this PR.
Done. Removed EvaluateGradient() from MeshField. It is now in MeshFieldLinear only. It also simplified to EvaluateGradient(element) instead of EvaluateGradient(element, barycentric).
geometry/proximity/mesh_field_linear.h, line 41 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW I'd happily swap
b_EQfor justb_Q(as I believe that's notation that I've used elsewhere (e.g.,surface_mesh.h, and most other places).
Done. PTAL.
There was a problem hiding this comment.
Strap in, here's a bunch of comments.
Reviewed 17 of 24 files at r1, 3 of 3 files at r2, 2 of 2 files at r3.
Reviewable status: 45 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy and @SeanCurtis-TRI)
a discussion (no related file):
I'm going to propose that we pull the contact_surface_doxygen.h out of this PR. I think there's a lot of word-smithing that could take place better in a google doc. I also want the chance to play with the figures. I feel that some points could be made better with illustrations instead of screen grabs of 3D scenes. By backing it out of this PR, we'll have more of an opportunity to refine it and it reduces the size of this PR significantly.
a discussion (no related file):
Make sure there's a test on mesh_intersection that the contact surface's gradients are reported in the world frame.
geometry/proximity/collision_doxygen.h, line 6 at r3 (raw file):
nit: I'd recommend the following:
Penetration queries report all pairs of intersecting geometries and differ in how they characterize the penetration. They include:
geometry/proximity/collision_doxygen.h, line 10 at r3 (raw file):
- report penetrations as point pairs - @ref module_contact_surface "report penetrations as contact surfaces" - check whether there is any collision at all
nit s/collision/penetration
geometry/proximity/contact_surface_doxygen.h, line 8 at r3 (raw file):
A contact surface is the surface of pressure equilibrium between two compliant objects in the hydroelastic contact model.
BTW If you mean rigid objects to be "compliant objects" than this is strictly true. But the language suggests something else. And what it suggests (soft) isn't strictly true for the code.
geometry/proximity/contact_surface_doxygen.h, line 10 at r3 (raw file):
compliant objects in the hydroelastic contact model. Rigid objects do not have pressure fields, so we construct the contact
BTW Here's an alternative way to say this which is consistent with the "equilibrium" surface idea.
A rigid object's pressure field goes from 0 pressure to infinite pressure instantaneously -- therefore the equilibrium surface must lie on the rigid object's surface.
It's also the surface as the relative "stiffness" of one of the surfaces goes to infinity.
What you write here depends on what level of abstraction you're targeting -- is this summary about the concept or about the implementation.
Although, I admit there's some sloppiness in that analogy.
geometry/proximity/contact_surface_doxygen.h, line 17 at r3 (raw file):
<h2> Thin rigid objects </h2> When a rigid object N deeply penetrates into a soft object M, their contact
BTW N = rigid and M = soft is less intuitive than rigid object R and soft object S.
geometry/proximity/contact_surface_doxygen.h, line 18 at r3 (raw file):
nit: A legalistic point: the contact surface with "undesirable parts" is only according to a strict definition. You should couch it as such. We don't want to suggest we're building contact surfaces with undesirable parts. :) I.e.,
If we define the contact surface strictly as the intersection of rigid surface mesh with soft volume mesh, the resulting contact surface would have undesirable artifacts.
geometry/proximity/contact_surface_doxygen.h, line 19 at r3 (raw file):
When a rigid object N deeply penetrates into a soft object M, their contact surface may have undesirable parts. For example, the following picture shows a rigid thin plate in transparent grey deeply penetrates into a soft ball
nit: s/penetrates/penetrating
geometry/proximity/contact_surface_doxygen.h, line 20 at r3 (raw file):
nit: "a thing shown in a color" is more directly expressed as "a color thing".
Image shows a grey rigid thin plane deeply penetrating into a pink soft ball.
(The transparency is probably sufficiently obvious.)
geometry/proximity/contact_surface_doxygen.h, line 20 at r3 (raw file):
surface may have undesirable parts. For example, the following picture shows a rigid thin plate in transparent grey deeply penetrates into a soft ball shown in transparent pink. The picture shows two surface patches of the
nit "...the picture shows the two surface..."
NOTE: I stopped reading at this point for the reason given above.
geometry/proximity/mesh_field_linear.h, line 23 at r2 (raw file):
Previously, DamrongGuoy (Damrong Guoy) wrote…
Done.
Yep. That does the summary job. Thanks
geometry/proximity/mesh_field_linear.h, line 64 at r2 (raw file):
Previously, DamrongGuoy (Damrong Guoy) wrote…
∇uᵉ depends on the values of u at vertices. ∇bᵢ does not depend on values of u.
oops. I msiread that. :) Withdrawn
geometry/proximity/mesh_field_linear.h, line 121 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: This shouldn't be public.
Per our f2f -- this proved to be more complex than at first glance.
Originally, a limiting factor for making this private is that we have mutable_values() which allows someone to change the Fi values but leaving the \nablaFi values unchanged -- leading to an inconsistent mesh field.
We also know that when a contact surface re-expresses itself (if it needs to switch ids) it goes from being expressed in one frame to another and these gradients would likewise have to be re-expressed. Otherwise, the mesh isn't consistent with the mesh.
Looking in mesh_intersection.h (line 801) we noted some additional difficulties. Here's our conclusion (as I see it).
- Delete the mutable_values(); it is only used once in
contact_surface_test.cc-- we can swap that use for something equivlant to achieve the same end. - Add a method to
MeshFieldLinear-- i likeReExpress(RigidTransform), but if you wanted to make something more literaly likeTransformGradients()to match theMesh::TransformVertices()I could accept that symmetry. - In the location indicated in
mesh_intersection.h, add the invocation to re-express the mesh field. - Make the gradient calculation private and call it in the constructor (this latter is already done).
- Add a TODO somewhere indicating that we want to compute the mesh and field with the quantities expressed in the target frame by construction so that we can delete those two transforming methods.
geometry/proximity/mesh_field_linear.h, line 28 at r3 (raw file):
BTW Related to a note below, perhaps this should be:
This class approximate the field
fasf*, by storing ...
Feel free to swap f* with whatever symbol you like indicating the function that this mesh field actually is. If you want to make it fᵉ I think you'll have the challenge to justify the e as mnemonic, but I wouldn't dismiss it out of hand.
geometry/proximity/mesh_field_linear.h, line 31 at r3 (raw file):
BTW I'd suggest one last bit to tie this all together:
Each element is defined by a tuple of index values
(i0, i1, i2)for tris,(i0, i1, i2, i3)for tets, such that for a given element E: V(i0) is the position of the zeroth vertex referenced by the element and F(i0) is the value offat that same position.
Or some such thing.
geometry/proximity/mesh_field_linear.h, line 44 at r3 (raw file):
BTW I'd suggest injecting a sentence:
simplicial element E. We indicate the barycentric coordinate representation of a point Q on an element E as b(Q) -- omitting the E for typographical convenience. The coefficients...
Or some such thing.
geometry/proximity/mesh_field_linear.h, line 46 at r3 (raw file):
to identify a point Q that lies in the simplicial element E. The coefficients b₀, b₁, b₂ (, b₃ for tetrahedron) are the weights of vertices V[i₀], V[i₁], V[i₂] (, V[i₃] for tetrahedron) of element E. Each bᵢ is positive or zero, and
BTW You might consider throwing out "Each bᵢ is positive or zero, and their summation is 1." because it's well defined just five lines higher in the definition of barycentric coordinates.
geometry/proximity/mesh_field_linear.h, line 49 at r3 (raw file):
their summation is 1. A point Q can be expressed as: Q = b₀(Q)V[i₀] + b₁(Q)V[i₁] + b₂(Q)V[i₂] for triangle,
BTW We probably want to introduce the symbols b_0(Q), i_0, etc. a bit more formally. The definition of b(Q) can simply be introduced in the discussion of barycentric coordinates. And the i_0 notation can in the summary -- a feature of the mesh.
geometry/proximity/mesh_field_linear.h, line 57 at r3 (raw file):
field f at Q is: fᵉ(Q) = b₀(Q)F[i₀] + b₁(Q)F[i₁] + b₂(Q)F[i₂] for triangle,
BTW
NOTE: After writing this note, I went up and modified/added notes higher in this file. They represent the implications of the position I'm arguing here. So, if this seems somewhat redundant in ways, that's the reason.
I'm not a fan of the super-script e on f. It doesn't have meaning for me. Given function f, this class represents f*, according to my sketchy notes, or L(f) based on some past emails. Evaluating that function at a point Q is just that: f*(Q) or L(f(Q)). Mathematically, the element on which Q lies is a hidden detail. It is, of course, very relevant to the computation (and is part of the API for reasons of optimization). But as this documentation is about the math, I don't think that API implementation detail should appear.
If there were, in fact, any e appearing in this math at all, it would be on the definition of the barycentric coordinate and the selection of the face that is providing the function indices. And I like the idea of omitting that in the mathematical expression for simplicity relying on the sentence leading up to the math to provide that context.
geometry/proximity/mesh_field_linear.h, line 60 at r3 (raw file):
= b₀(Q)F[i₀] + b₁(Q)F[i₁] + b₂(Q)F[i₂] + b₃(Q)F[i₃] for tetrahedron, where F[i₀], F[i₁], F[i₂] (, F[i₃] for tetrahedron) are the field values at
BTW If you take the earlier decision in the summary where the i0 notation is introduced, you don't need to re-define those symbols here.
geometry/proximity/mesh_field_linear.h, line 65 at r3 (raw file):
nit: This sentence
The gradient ∇bᵢ is constant on E.
is duplicated after the math definition. Can we streamline this to say it only once in such close proximity?
geometry/proximity/mesh_intersection.h, line 519 at r3 (raw file):
● ↗ ↑ ↖ ● ↗ pressure gradient ∇p_M in M ● ↗ ⇧ ↑ ⇧ ↖ ● ⇧ face normal f_N of N ● ┌───────┐ ●
Probably want to apply the "suitable"/"unsuitable" lines on this one as well. Although, in this case the point would be what was "deemed" unsuitable isn't really.
geometry/proximity/mesh_intersection.h, line 547 at r3 (raw file):
*/ template <typename T> bool IsFaceNormalAlongPressureGradient(
BTW The name of this method is very literal. It describes the operation instead of the semantics. You might consider renaming it to something that describes more what is being accomplished. Perhaps:
MeshNormalCompatibleWithGradient() or
TriangleNeedsCulling() (this would require reversing the semantics but has the advantage of reversing the the negated boolean expressions in the invocations below).
geometry/proximity/mesh_intersection.h, line 563 at r3 (raw file):
// We pick 5π/8 empirically. const double kAngleThreshold = 5. * M_PI / 8.; const T kDotProductThreshold = std::cos(kAngleThreshold);
BTW Because this isn't a constexpr (and can't be because std::cos isn't declared as such), you might consider pre-computing this somewhere once and re-using it, or pre-computing offline and simply setting the value to the literal value.
However, if you feel that is a pre-optimization and want to just leave it, as is, then how about a TODO to consider removing it?
geometry/proximity/proximity_doxygen.h, line 8 at r3 (raw file):
Proximity queries span a range of types, including: - @ref module_collision_queries "collision"
nit: s/collision/penetration (in both the doxygen tag name and the string.
If you look at the group, you'll see "Penetration" everywhere. The function names are called that, the return types are called that. Let's go ahead and be consistent.
geometry/proximity/proximity_doxygen.h, line 20 at r3 (raw file):
@{ @defgroup module_collision_queries Collision Queries
nit: s/Collision/Penetration
geometry/proximity/surface_mesh.h, line 106 at r3 (raw file):
}; // Forward declaration of SurfaceMeshTester<T>. SurfaceMesh<T> will
nit: It would be nice to dispose of the forward declarations here and in the volume_mesh.h.
geometry/proximity/surface_mesh.h, line 133 at r3 (raw file):
A surface mesh has the intrinsic dimension 2. It is embedded in 3-d space. */ static constexpr int kDim = 2;
kDim largely seems pointless now. Most places in geometry/ where we use it, is in a kDim + 1 kind of way. So, we might as well replace kDim with kVertexPerElement and keep things simpler.
geometry/proximity/surface_mesh.h, line 134 at r3 (raw file):
*/ static constexpr int kDim = 2; static constexpr int kVertexPerElement = 3;
nit: This is the only type trait lacking doxygen.
geometry/proximity/surface_mesh.h, line 426 at r3 (raw file):
const std::array<FieldValue, 3>& field_value, SurfaceFaceIndex f) const { Vector3<FieldValue> gradu_M = field_value[0] * CalcGradBarycentric(f, 0); for (int i = 1; i < 3; ++i) {
BTW this function would be slightly shorter if you unrolled the for loop by hand. E.g.,
Vector3<FieldValue> gradu_M = field_value[0] * CalcGradBarycentric(f, 0);
gradu_M += field_value[1] * CalcGradBarycentric(f, 1);
gradu_M += field_value[2] * CalcGradBarycentric(f, 2);geometry/proximity/surface_mesh.h, line 503 at r3 (raw file):
Vector3<T> SurfaceMesh<T>::CalcGradBarycentric(SurfaceFaceIndex f, int i) const { DRAKE_DEMAND(0 <= i && i < 3);
This function would definitely benefit from some implementation documentation explaining the math.
geometry/proximity/surface_mesh.h, line 504 at r3 (raw file):
int i) const { DRAKE_DEMAND(0 <= i && i < 3); const Vector3<T>& V = vertices_[faces_[f].vertex(i)].r_MV();
nit: How about some monogram notation?
A-->p_AMB-->p_BM
geometry/proximity/volume_mesh.h, line 130 at r3 (raw file):
static constexpr int kDim = 3; static constexpr int kVertexPerElement = 4;
nit: Same note as for surface mesh.
geometry/proximity/volume_mesh.h, line 256 at r3 (raw file):
@returns `true` if the given mesh is equal. */ bool Equal(const VolumeMesh<T>& mesh) const {
BTW this isn't your fault and it isn't new to this PR, but this Equal() method seems to have a defect (same in SurfaceMesh).
If I do mesh.Equal(mesh);, it should return true trivially.
geometry/proximity/volume_mesh.h, line 311 at r3 (raw file):
int i) const { DRAKE_DEMAND(0 <= i && i < 4); const Vector3<T>& V = vertices_[elements_[e].vertex(i)].r_MV();
nit: Same notes as from the SurfaceMesh version of this function.
geometry/proximity/test/mesh_field_linear_test.cc, line 107 at r3 (raw file):
} GTEST_TEST(MeshFieldLinearTest, TestEvaluateGradient) {
BTW We should probably test this with VolumeMesh. Imagine this were mistakenly hard-coded with for loops 0 <= i < 3;. This test would pass but volume mesh would fail. So, let's do both.
geometry/proximity/test/mesh_intersection_test.cc, line 408 at r3 (raw file):
unique_ptr<VolumeMeshFieldLinear<T, T>> TrivialVolumeMeshField( const VolumeMesh<T>* volume_mesh) { // TODO(SeanCurtis-TRI): All the zeros and ones prevent meaningful recognition
nit: This TODO is still relevant; scaling the ones does not change the concern.
geometry/proximity/test/mesh_intersection_test.cc, line 678 at r3 (raw file):
You should add a comment along the lines of:
It is ok to use the trivial mesh and trivial mesh field in this test. The function under test asks for the gradient values and operates on it. It is not responsible for making sure that the gradient is computed correctly -- that is tested elsewhere.
geometry/proximity/test/mesh_intersection_test.cc, line 688 at r3 (raw file):
// We will set the pose of SurfaceMesh N in M's frame so that the triangle // Face_0 of N has its face normal vector make various angle with the
nit s/angle/angles
geometry/proximity/test/mesh_intersection_test.cc, line 693 at r3 (raw file):
double angle; // Angle between the face normal and the gradient. bool expect_result; // true when `angle` < threshold 5π/8. } test_data[6]{{0, true},
nit: You can omit the size (6) and let the compiler infer it from initialization.
geometry/proximity/test/mesh_intersection_test.cc, line 707 at r3 (raw file):
// the triangle and the tetrahedron intersect. const auto X_MN = RigidTransformd(RollPitchYawd(t.angle, 0, 0), Vector3d::Zero());
BTW This is my broken drum. :) If you rotated around some arbitrary vector, your test would be more robust. (As long as its not Mz, you're fine.)
geometry/proximity/test/surface_mesh_test.cc, line 334 at r3 (raw file):
const Vector3<T> expect_gradb2_M = Vector3<T>::UnitY() / 2.; const auto R_WM = X_WM.rotation();
nit const auto&
geometry/proximity/test/surface_mesh_test.cc, line 339 at r3 (raw file):
EXPECT_TRUE(CompareMatrices(R_WM * expect_gradb2_M, gradb2_W, 1e-14)); // Check the invariance ∑(∇bᵢ) = 0.
nit = 1.
geometry/proximity/test/surface_mesh_test.cc, line 344 at r3 (raw file):
} GTEST_TEST(SurfaceMeshTest, TestCalcGradBarycentricDouble) {
BTW This is what typed parameter tests are for. :) However, you don't have enough critical mass yet to require it.
geometry/proximity/test/volume_mesh_test.cc, line 103 at r3 (raw file):
TestVolumeMesh<T>(RigidTransform<T>(Vector3<T>(1., 1., 1.))); EXPECT_TRUE(mesh->Equal(*mesh_too)); EXPECT_FALSE(mesh->Equal(*mesh_translate));
nit: As for testing Equal, you've only tested one way for them to be unequal. There are lots of expected conditions that would invalidate equality; you might consider either TODO'ing those or DO'ing those. :)
geometry/proximity/test/volume_mesh_test.cc, line 262 at r3 (raw file):
EXPECT_TRUE(CompareMatrices(R_WM * expect_gradb3_M, gradb3_W, 1e-14)); // Check the invariance ∑(∇bᵢ) = 0.
nit = 1.
geometry/proximity/test/volume_mesh_test.cc, line 355 at r3 (raw file):
Vector3<T> gradp_M = volume_mesh_field->EvaluateGradient(e0, b); Vector3<T> expect_gradp_M = volume_mesh->CalcGradientVectorOfLinearField( std::array<T, 4>{p0, p1, p2, p3}, e0);
BTW In all of these tests, you might want to swap pi for fi or ui. p has "point-like" connotations and your documentation uses f. SO, might as well tie them together better.
There was a problem hiding this comment.
Reviewed 7 of 7 files at r4.
Reviewable status: 45 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy)
geometry/proximity/mesh_field_linear.h, line 125 at r4 (raw file):
Let's try this:
Evaluates the gradient in the domain of the element indicated by
e. The gradient is a vector in R3 expressed in frame M. For surface meshes, it will particularly lie parallel to the plane of the corresponding triangle.
Feel free to rephrase as you like.
There was a problem hiding this comment.
Reviewable status: 39 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @DamrongGuoy and @SeanCurtis-TRI)
a discussion (no related file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
I'm going to propose that we pull the contact_surface_doxygen.h out of this PR. I think there's a lot of word-smithing that could take place better in a google doc. I also want the chance to play with the figures. I feel that some points could be made better with illustrations instead of screen grabs of 3D scenes. By backing it out of this PR, we'll have more of an opportunity to refine it and it reduces the size of this PR significantly.
I moved contact_surface_doxygen.h and its images to Goodgle Drive.
https://drive.google.com/drive/folders/1mgsHdj0kBf8UoR7zfzDHfax2ebLvLLoR
geometry/proximity/mesh_field_linear.h, line 121 at r2 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
Per our f2f -- this proved to be more complex than at first glance.
Originally, a limiting factor for making this private is that we have
mutable_values()which allows someone to change the Fi values but leaving the \nablaFi values unchanged -- leading to an inconsistent mesh field.We also know that when a contact surface re-expresses itself (if it needs to switch ids) it goes from being expressed in one frame to another and these gradients would likewise have to be re-expressed. Otherwise, the mesh isn't consistent with the mesh.
Looking in mesh_intersection.h (line 801) we noted some additional difficulties. Here's our conclusion (as I see it).
- Delete the mutable_values(); it is only used once in
contact_surface_test.cc-- we can swap that use for something equivlant to achieve the same end.- Add a method to
MeshFieldLinear-- i likeReExpress(RigidTransform), but if you wanted to make something more literaly likeTransformGradients()to match theMesh::TransformVertices()I could accept that symmetry.- In the location indicated in
mesh_intersection.h, add the invocation to re-express the mesh field.- Make the gradient calculation private and call it in the constructor (this latter is already done).
- Add a TODO somewhere indicating that we want to compute the mesh and field with the quantities expressed in the target frame by construction so that we can delete those two transforming methods.
Done.
- Removed
mutable_values()and used slightly more code incontact_surface_test.cc. - Added
MeshFieldLinear::TransformGradients(). - In
mesh_intersection.h, added calls toTransformGradients()whenTransformVertices()were called. CalcGradientFieldis now private and still called by the contructor.- Added TODO about the target World frame in
mesh_intersection.hwhere we callTransformGradients()andTransformVertices().
geometry/proximity/mesh_intersection.h, line 519 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
Probably want to apply the "suitable"/"unsuitable" lines on this one as well. Although, in this case the point would be what was "deemed" unsuitable isn't really.
Done. I changed the word "unsuitable" to "incorrectly prohibited".
geometry/proximity/mesh_intersection.h, line 547 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW The name of this method is very literal. It describes the operation instead of the semantics. You might consider renaming it to something that describes more what is being accomplished. Perhaps:
MeshNormalCompatibleWithGradient()or
TriangleNeedsCulling()(this would require reversing the semantics but has the advantage of reversing the the negated boolean expressions in the invocations below).
I'd like to keep the name at least for now.
geometry/proximity/mesh_intersection.h, line 563 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW Because this isn't a
constexpr(and can't be becausestd::cosisn't declared as such), you might consider pre-computing this somewhere once and re-using it, or pre-computing offline and simply setting the value to the literal value.However, if you feel that is a pre-optimization and want to just leave it, as is, then how about a TODO to consider removing it?
Done. I added a TODO. It is also possible that users might need it as an argument to the function.
geometry/proximity/proximity_doxygen.h, line 8 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: s/collision/penetration (in both the doxygen tag name and the string.
If you look at the group, you'll see "Penetration" everywhere. The function names are called that, the return types are called that. Let's go ahead and be consistent.
I got "Collision" from query_object.h:
$ grep Collision geometry/query_object.h
@name Collision Queries
std::vector<SortedPair<GeometryId>> FindCollisionCandidates() const;
bool HasCollisions() const;
I also see "Penetration" in query_object.h:
$ grep Penetration geometry/query_object.h
double-valued scalar -- see ComputePointPairPenetration()). In other cases,
PenetrationAsPointPair), providing some measure of the penetration "depth" of
std::vector<PenetrationAsPointPair<double>> ComputePointPairPenetration()
I changed to "penetration".
geometry/proximity/proximity_doxygen.h, line 20 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: s/Collision/Penetration
Done.
geometry/proximity/collision_doxygen.h, line 6 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: I'd recommend the following:
Penetration queries report all pairs of intersecting geometries and differ in how they characterize the penetration. They include:
Done.
geometry/proximity/collision_doxygen.h, line 10 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit s/collision/penetration
Done.
There was a problem hiding this comment.
Reviewable status: 34 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @DamrongGuoy and @SeanCurtis-TRI)
geometry/proximity/mesh_field_linear.h, line 28 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW Related to a note below, perhaps this should be:
This class approximate the field
fasf*, by storing ...Feel free to swap
f*with whatever symbol you like indicating the function that this mesh field actually is. If you want to make itfᵉI think you'll have the challenge to justify theeas mnemonic, but I wouldn't dismiss it out of hand.
I use fᵉ for the part of the Finite Element Approximation on element E. It is somewhat standard in Finite Element literature. I never see f*.
geometry/proximity/mesh_field_linear.h, line 31 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW I'd suggest one last bit to tie this all together:
Each element is defined by a tuple of index values
(i0, i1, i2)for tris,(i0, i1, i2, i3)for tets, such that for a given element E: V(i0) is the position of the zeroth vertex referenced by the element and F(i0) is the value offat that same position.Or some such thing.
I'd like to think of vertices and elements in abstract sense without too much into indexing.
geometry/proximity/mesh_field_linear.h, line 46 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW You might consider throwing out "Each bᵢ is positive or zero, and their summation is 1." because it's well defined just five lines higher in the definition of barycentric coordinates.
I'd like to keep both math equations and english words.
There was a problem hiding this comment.
Reviewable status: 34 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @DamrongGuoy and @SeanCurtis-TRI)
geometry/proximity/mesh_field_linear.h, line 57 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW
NOTE: After writing this note, I went up and modified/added notes higher in this file. They represent the implications of the position I'm arguing here. So, if this seems somewhat redundant in ways, that's the reason.
I'm not a fan of the super-script
eonf. It doesn't have meaning for me. Given functionf, this class representsf*, according to my sketchy notes, orL(f)based on some past emails. Evaluating that function at a point Q is just that:f*(Q)orL(f(Q)). Mathematically, the element on which Q lies is a hidden detail. It is, of course, very relevant to the computation (and is part of the API for reasons of optimization). But as this documentation is about the math, I don't think that API implementation detail should appear.If there were, in fact, any
eappearing in this math at all, it would be on the definition of the barycentric coordinate and the selection of the face that is providing the function indices. And I like the idea of omitting that in the mathematical expression for simplicity relying on the sentence leading up to the math to provide that context.
Now I think I see the miscommunication. I'm not trying to use fᵉ for f*. In Finite Element literature, fᵉ is defined only on the element E, i.e., fᵉ : E -> R. I think your f* is defined on the whole mesh M, i.e., f* : M -> R. I have not tried to use any symbol for the function M->R. I'll have to look up FEM book for that symbol, perhaps f̃.
I try to focus the discussion only on one element since it's already long enough.
Every time I mention fᵉ, I start with "On the element E". I'll try to make it clearer.
DamrongGuoy
force-pushed
the
field_gradient
branch
from
d3e6a10
to
1f724e5
Compare
There was a problem hiding this comment.
Reviewed 10 of 17 files at r5.
Reviewable status: 34 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy)
geometry/proximity/mesh_field_linear.h, line 28 at r3 (raw file):
While none of this documentation alludes to finite elements, I can accept that. However, you still need to introduce this symbol up here:
This class approximates the field f as fᵉ...
geometry/proximity/mesh_field_linear.h, line 31 at r3 (raw file):
Previously, DamrongGuoy (Damrong Guoy) wrote…
I'd like to think of vertices and elements in abstract sense without too much into indexing.
That's fine. You get to think about it any way you like. :) However, this is notation you use below that isn't defined up here. I tend to balk at ad hoc notation that simply appears. If you want to use it, it's got to be defined.
geometry/proximity/mesh_field_linear.h, line 57 at r3 (raw file):
I really like this point.
But there's still a subtle issue I have with it. I think the symbol f is still overloaded. It lends itself to the interpretation of evaluating f in the domain of E. But that's obviously not the case. We've encountered problems in the past with using the symbol for the original function and the linear-approximation of the function interchangeably and that led to a mountain of confusion and error. That is the issue I'm trying to avoid.
However, if, as I've previously indicated, you introduced this symbol up above, I can accept it. E.g.,
We approximate the field
fwith a set of piecewise linear functions -- onefᵉfor each elementEsuch that the aggregate effect is a continuous piecewise linear approximation of f.
(Or some such thing. I struggled with the second clause to efficiently communicate that the linear functions are defined such that they "line up" at common domain boundaries. I'll leave that up to you.)
geometry/proximity/mesh_field_linear.h, line 138 at r5 (raw file):
const math::RigidTransform<typename MeshType::ScalarType>& X_NM) { for (auto& grad : gradients_) { // Alias `grad` as expressed in frame M.
BTW This seems like a great deal of counter-intuitive work. In this case, I'd just go ahead with:
for (auto& grad : gradients_) {
grad = X_NM.rotation() * grad;
}I think in this case, the amount of work is so clear and obvious that you don't risk anything notationally, and the code will be much more obvious in cause and effect. @sherm1 what do you think?
geometry/proximity/proximity_doxygen.h, line 8 at r3 (raw file):
Previously, DamrongGuoy (Damrong Guoy) wrote…
I got "Collision" from query_object.h:
$ grep Collision geometry/query_object.h @name Collision Queries std::vector<SortedPair<GeometryId>> FindCollisionCandidates() const; bool HasCollisions() const;I also see "Penetration" in query_object.h:
$ grep Penetration geometry/query_object.h double-valued scalar -- see ComputePointPairPenetration()). In other cases, PenetrationAsPointPair), providing some measure of the penetration "depth" of std::vector<PenetrationAsPointPair<double>> ComputePointPairPenetration()I changed to "penetration".
Thanks. :)
geometry/proximity/surface_mesh.h, line 511 at r5 (raw file):
const Vector3<T> ehat = (B - A).normalized(); const Vector3<T> AV = V - A; const Vector3<T> nhat = (AV - (AV.dot(ehat)) * ehat).normalized();
nit: You'll want to document why this function is ok if the point V is co-linear with the AB edge.
geometry/proximity/test/contact_surface_test.cc, line 48 at r5 (raw file):
using Eigen::Vector3d; using std::make_unique; using std::move;
nit: By adding these aliases, you've left the file in a somewhat incoherent state. We still have lots of usages of std::make_unique and std::move. If you add the aliases, you should unify the file. Perhaps better for a different)PR. IN this case, it might be better to simply use std:: in the few lines of code you've added below.
geometry/proximity/test/mesh_field_linear_test.cc, line 121 at r5 (raw file):
"e", std::move(e_values), mesh.get()); Vector3<double> gradient = mesh_field->EvaluateGradient(SurfaceFaceIndex(0));
nit: Now that you've aliased Vector3d could you swap this (and the next) declarations to use that? In this case, they are the only two outliers affected by the new alias. Worth doing now.
There was a problem hiding this comment.
Reviewed 1 of 24 files at r1, 1 of 7 files at r4, 2 of 17 files at r5.
Reviewable status: 45 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy)
a discussion (no related file):
@DamrongGuoy, I believe you might want to take a look at "Hughes, 2000. The Finite Element Method", chapter 3.9 (I believe you bought this book when you asked me about FEM). Don't get scared, that chapter is only a couple pages long. It shows how FEM people compute the gradients of the shape functions (your CalcGradBarycentric() method). Your implementation is great! I love the geometric way of doing this (vs. the FEM algebraic way). However I wonder, the fact we choose a particular edge (in your case you chose AV, why not BV?) as a reference, would it have implications regarding accuracy? especially for thing triangles as those in hydroelastic that get formed as a trinagle first shows up and grows from zero size.
geometry/proximity/mesh_field.h, line 101 at r5 (raw file):
extern template class MeshField<double, SurfaceMesh<double>>; extern template class MeshField<AutoDiffXd, SurfaceMesh<AutoDiffXd>>;
nit, why not DRAKE_DECLARE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS?
geometry/proximity/mesh_field.cc, line 7 at r5 (raw file):
template class MeshField<double, SurfaceMesh<double>>; template class MeshField<AutoDiffXd, SurfaceMesh<AutoDiffXd>>;
nit, why not DRAKE_DEFINE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS()
geometry/proximity/mesh_field_linear.h, line 136 at r5 (raw file):
*/ void TransformGradients( const math::RigidTransform<typename MeshType::ScalarType>& X_NM) {
sry if I do not understand the use case. Isn't it dangerous we can mutate the field's gradients without somehow making sure it still is consistent with the original mesh M it makes reference to? maybe at least an "(advanced)" or "@warning" is in order?
geometry/proximity/mesh_intersection.h, line 568 at r5 (raw file):
// We pick 5π/8 empirically. const double kAngleThreshold = 5. * M_PI / 8.;
why not constexpr? then you could resolve the TODO above.
geometry/proximity/surface_mesh.h, line 507 at r5 (raw file):
DRAKE_DEMAND(0 <= i && i < 3); const Vector3<T>& V = vertices_[faces_[f].vertex(i)].r_MV(); const Vector3<T>& A = vertices_[faces_[f].vertex((i + 1) % 3)].r_MV();
nit, use monogram. i.e. A --> p_MA and B --> p_MB
geometry/proximity/surface_mesh.h, line 508 at r5 (raw file):
const Vector3<T>& V = vertices_[faces_[f].vertex(i)].r_MV(); const Vector3<T>& A = vertices_[faces_[f].vertex((i + 1) % 3)].r_MV(); const Vector3<T>& B = vertices_[faces_[f].vertex((i + 2) % 3)].r_MV();
nit, why not p_MV0, p_MV1,p_MV2 instead of V, A and B?
geometry/proximity/surface_mesh.h, line 509 at r5 (raw file):
const Vector3<T>& A = vertices_[faces_[f].vertex((i + 1) % 3)].r_MV(); const Vector3<T>& B = vertices_[faces_[f].vertex((i + 2) % 3)].r_MV(); const Vector3<T> ehat = (B - A).normalized();
nit, use monogram ehat --> ehat_M
geometry/proximity/surface_mesh.h, line 510 at r5 (raw file):
const Vector3<T>& B = vertices_[faces_[f].vertex((i + 2) % 3)].r_MV(); const Vector3<T> ehat = (B - A).normalized(); const Vector3<T> AV = V - A;
nit, maybe p_V0V1?
geometry/proximity/surface_mesh.h, line 511 at r5 (raw file):
const Vector3<T> ehat = (B - A).normalized(); const Vector3<T> AV = V - A; const Vector3<T> nhat = (AV - (AV.dot(ehat)) * ehat).normalized();
I believe you can get away with a single norm (less expensive):
e = B - A;
AV = V - A;
n = AV - AV.dot(e) * e / e.squaredNorm();
grad_Ni = n / n.norm();
return grad_Ni;where I use the notation Ni for the i-th shape function (as it is traditionally done for finite elements for instance, though not unique to finite elements but to the geometry of it).
geometry/proximity/surface_mesh.h, line 512 at r5 (raw file):
const Vector3<T> AV = V - A; const Vector3<T> nhat = (AV - (AV.dot(ehat)) * ehat).normalized(); // nhat.dot(AV) = nhat.dot(BV)
is this true? consider a right triangle with BV perpendicular to AB (i.e. BV aligned with nhat).
DamrongGuoy
force-pushed
the
field_gradient
branch
2 times, most recently
from
039f6e1
to
d004a7b
Compare
There was a problem hiding this comment.
Reviewable status: 41 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @amcastro-tri, @DamrongGuoy, @SeanCurtis-TRI, and @warning)
a discussion (no related file):
Previously, amcastro-tri (Alejandro Castro) wrote…
@DamrongGuoy, I believe you might want to take a look at "Hughes, 2000. The Finite Element Method", chapter 3.9 (I believe you bought this book when you asked me about FEM). Don't get scared, that chapter is only a couple pages long. It shows how FEM people compute the gradients of the shape functions (your
CalcGradBarycentric()method). Your implementation is great! I love the geometric way of doing this (vs. the FEM algebraic way). However I wonder, the fact we choose a particular edge (in your case you chose AV, why not BV?) as a reference, would it have implications regarding accuracy? especially for thing triangles as those in hydroelastic that get formed as a trinagle first shows up and grows from zero size.
Thank you for pointing out to Hughes's book. I have just read it.
Indeed I previously formulated the solution in terms of matrix inversion similar to Hughes's book, which is more general and applies to many kinds of elements (simplicial, cuboid,...). However, I found that the specific geometric solution that I used was easier to understand, though less general.
Both methods suffer from thin triangles. The matrix inversion would have a poor-conditioned matrix. The geometric formula also suffer.
Similar to the older PR about setting face normals of contact surfaces (SurfaceMesh::CalcAreasNormalsAndCentroid()), I put a TODO to watch out for almost-zero area triangles in SurfaceMesh::CalcGradBarycentric(). Similar to the TODO to set face normals on contact polygons from face normals of rigid triangles, we can set the pressure gradient along a contact polygon from the projection of the soft tetrahedron's pressure gradient onto the rigid triangle.
If it sounds reasonable, please unblock.
geometry/proximity/mesh_field.h, line 101 at r5 (raw file):
Previously, amcastro-tri (Alejandro Castro) wrote…
nit, why not
DRAKE_DECLARE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS?
Because MeshField template takes two parameters: FieldValue and MeshType. The DRAKE_... macro only allows one parameter.
geometry/proximity/mesh_field.cc, line 7 at r5 (raw file):
Previously, amcastro-tri (Alejandro Castro) wrote…
nit, why not
DRAKE_DEFINE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS()
Because MeshField template takes two parameters: FieldValue and MeshType. The DRAKE_... macro only allows one parameter.
geometry/proximity/mesh_field_linear.h, line 28 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
While none of this documentation alludes to finite elements, I can accept that. However, you still need to introduce this symbol up here:
This class approximates the field f as fᵉ...
Now I simplify the discussion by saying from the beginning that f itself is a piecewise linear field. It eliminates the need for fᵉ. Now I can see that fᵉ can confuse readers.
geometry/proximity/mesh_field_linear.h, line 44 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW I'd suggest injecting a sentence:
simplicial element E. We indicate the barycentric coordinate representation of a point Q on an element E as b(Q) -- omitting the E for typographical convenience. The coefficients...
Or some such thing.
Done six lines below. That's a good addition.
geometry/proximity/mesh_field_linear.h, line 49 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW We probably want to introduce the symbols
b_0(Q),i_0, etc. a bit more formally. The definition ofb(Q)can simply be introduced in the discussion of barycentric coordinates. And thei_0notation can in the summary -- a feature of the mesh.
Now I mention that i is local index, not global index. (I found that "index of index" i₀, i₁, i₂ is less clear to me.)
Local indexing is a common practice in Finite Element Analysis. Sometimes they use i for local index and a for global index: vertex Vᵢ is the i-th vertex of the element; vertex Vₐ is the a-th vertex of the entire mesh.
Switching back to local index, I can use the summation ∑ and make the equations easier to read.
Furthermore the code (Evaluate(), CalcGradientVector()) also use local index i ∈ [0,1,...,kVertexPerElement).
geometry/proximity/mesh_field_linear.h, line 125 at r4 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
Let's try this:
Evaluates the gradient in the domain of the element indicated by
e. The gradient is a vector in R3 expressed in frame M. For surface meshes, it will particularly lie parallel to the plane of the corresponding triangle.Feel free to rephrase as you like.
Done. Thanks. I like its compactness.
geometry/proximity/mesh_field_linear.h, line 136 at r5 (raw file):
Previously, amcastro-tri (Alejandro Castro) wrote…
sry if I do not understand the use case. Isn't it dangerous we can mutate the field's gradients without somehow making sure it still is consistent with the original mesh M it makes reference to? maybe at least an "(advanced)" or "@warning" is in order?
Done adding @warning. Indeed, we have operations like SurfaceMesh::TransformVertices() to change the frame of ContactSurface. In the future, we might have VolumeMesh::TransformVertices() as well when we do soft-soft contacts. Both the mesh and the mesh field have to sync.
geometry/proximity/mesh_field_linear.h, line 138 at r5 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW This seems like a great deal of counter-intuitive work. In this case, I'd just go ahead with:
for (auto& grad : gradients_) { grad = X_NM.rotation() * grad; }I think in this case, the amount of work is so clear and obvious that you don't risk anything notationally, and the code will be much more obvious in cause and effect. @sherm1 what do you think?
Done.
geometry/proximity/mesh_intersection.h, line 568 at r5 (raw file):
Previously, amcastro-tri (Alejandro Castro) wrote…
why not
constexpr? then you could resolve the TODO above.
The short answer is std::cos cannot be a constexpr.
Gcc compiled this one ok:
constexpr double kAngleThreshold = 5. * M_PI / 8.;
constexpr double kDotProductThreshold = std::cos(kAngleThreshold);but CLion drew a red line under kDotProductThreshold with non-constexpr function cos cannot be used in a constant expresssion.
CLang gave this compile error error: constexpr variable kDotProductThreshold must be initialized by a constant expression.
I added the explanation in TODO that std::cos() cannot be a constant expression.
|
geometry/proximity/mesh_field_linear.h, line 136 at r5 (raw file): Previously, DamrongGuoy (Damrong Guoy) wrote…
I have just realized that |
There was a problem hiding this comment.
Reviewable status: 40 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @amcastro-tri, @DamrongGuoy, @SeanCurtis-TRI, and @warning)
a discussion (no related file):
By the way, after I rebased on the new master today, I found that #12748 has re-order entries in BUILD.bazel alphabetically.
Please do not look at the diff between r5 and r6 of BUILD.bazel, since it will look confusing. Instead, please look at the diff between r0 and the latest.
There was a problem hiding this comment.
Reviewable status: 39 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @amcastro-tri, @DamrongGuoy, and @SeanCurtis-TRI)
a discussion (no related file):
The matrix inversion would have a poor-conditioned matrix
This is difficult for me to see, do you have any insights? Definitely the inverse goes to infinity given the area goes to zero, though not clear that the condition number goes up.
geometry/proximity/mesh_intersection.h, line 568 at r5 (raw file):
Previously, DamrongGuoy (Damrong Guoy) wrote…
The short answer is
std::coscannot be aconstexpr.Gcc compiled this one ok:
constexpr double kAngleThreshold = 5. * M_PI / 8.; constexpr double kDotProductThreshold = std::cos(kAngleThreshold);but CLion drew a red line under
kDotProductThresholdwithnon-constexpr function cos cannot be used in a constant expresssion.CLang gave this compile error
error: constexpr variable kDotProductThreshold must be initialized by a constant expression.I added the explanation in TODO that
std::cos()cannot be a constant expression.
I don't understand why. @SeanCurtis-TRI, is this a bug in the compiler? I thought this was a very well supported c++11 feature?
There was a problem hiding this comment.
Reviewable status: 39 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @amcastro-tri, @DamrongGuoy, and @SeanCurtis-TRI)
geometry/proximity/mesh_intersection.h, line 568 at r5 (raw file):
Previously, amcastro-tri (Alejandro Castro) wrote…
I don't understand why. @SeanCurtis-TRI, is this a bug in the compiler? I thought this was a very well supported c++11 feature?
Anything declared to be constexpr must be initialized only with things that are likewise constexpr. While functions can be declared as constexpr, std::cos() is not one of those (https://en.cppreference.com/w/cpp/numeric/math/cos). Therefore, while the angle threshold can be constexpr the cut-off value kDotProductThreshold cannot be. Hence my concern that in this leaf function we redundantly evaluate the cosine of a constant over and over.
There was a problem hiding this comment.
Reviewable status: 39 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @amcastro-tri, @DamrongGuoy, and @SeanCurtis-TRI)
geometry/proximity/mesh_intersection.h, line 568 at r5 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
Anything declared to be
constexprmust be initialized only with things that are likewiseconstexpr. While functions can be declared asconstexpr,std::cos()is not one of those (https://en.cppreference.com/w/cpp/numeric/math/cos). Therefore, while the angle threshold can beconstexprthe cut-off valuekDotProductThresholdcannot be. Hence my concern that in this leaf function we redundantly evaluate the cosine of a constant over and over.
Two quick items:
static const double kDotProductThreshold = std::cos(kAngleThreshold);- Declaring it
staticmeans it will only get evaluated once. Solves my concern about redundant work. - The threshold should be of type
T, it should be of typedouble. It's just a constant and not a parameter that in any way depends on the scalar type.
There was a problem hiding this comment.
Reviewable status: 38 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @amcastro-tri, @DamrongGuoy, and @SeanCurtis-TRI)
geometry/proximity/mesh_intersection.h, line 568 at r5 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
Two quick items:
static const double kDotProductThreshold = std::cos(kAngleThreshold);
- Declaring it
staticmeans it will only get evaluated once. Solves my concern about redundant work.- The threshold should be of type
T, it should be of typedouble. It's just a constant and not a parameter that in any way depends on the scalar type.
TIL, I thought the math library supported constexprs. That's a shame.
BTW, I do like the static solution.
There was a problem hiding this comment.
Reviewable status: 38 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @amcastro-tri, @DamrongGuoy, and @SeanCurtis-TRI)
geometry/proximity/mesh_intersection.h, line 568 at r5 (raw file):
BTW My item 2 up above is missing a critical word that changes the entire meaning of the sentence:
The threshold should not be of type
T.
Just like me to drop the one word that matters most. :)
There was a problem hiding this comment.
Reviewed 4 of 4 files at r6.
Reviewable status: 38 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy)
geometry/proximity/mesh_field_linear.h, line 28 at r3 (raw file):
Previously, DamrongGuoy (Damrong Guoy) wrote…
Now I simplify the discussion by saying from the beginning that
fitself is a piecewise linear field. It eliminates the need forfᵉ. Now I can see thatfᵉcan confuse readers.
I love this. You are right to not worry about the fact that it's a linear approximation of anything. It's simpler to say, hey, it's a function and it's piecewise linear without commenting on why it exists. Two thumbs way up.
geometry/proximity/mesh_field_linear.h, line 49 at r3 (raw file):
Previously, DamrongGuoy (Damrong Guoy) wrote…
Now I mention that
iis local index, not global index. (I found that "index of index"i₀, i₁, i₂is less clear to me.)Local indexing is a common practice in Finite Element Analysis. Sometimes they use
ifor local index andafor global index: vertex Vᵢ is the i-th vertex of the element; vertex Vₐ is the a-th vertex of the entire mesh.Switching back to local index, I can use the summation
∑and make the equations easier to read.Furthermore the code (Evaluate(), CalcGradientVector()) also use local index i ∈ [0,1,...,kVertexPerElement).
Seems like a good compromise.
geometry/proximity/mesh_intersection.h, line 565 at r6 (raw file):
// std::cos(kAngleThreshold) once outside this function, perhaps as a // global variable. It is also possible that users might need this threshold // as an argument to the function. std::cos() cannot be a constant
BTW For fear of you putting you in the middle of a tug of war, I'd suggest this newly added sentence doesn't provide much value. I have several reasons, but they're not worth typing. Let me know if you want to hear the reasons.
There was a problem hiding this comment.
Reviewable status: 38 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy)
a discussion (no related file):
Previously, amcastro-tri (Alejandro Castro) wrote…
The matrix inversion would have a poor-conditioned matrix
This is difficult for me to see, do you have any insights? Definitely the inverse goes to infinity given the area goes to zero, though not clear that the condition number goes up.
Re the matrix condition number. Looking at Eq. 3.9.9 of Hugues book it does seem like the Jacobian, times its determinants (which equals 2 times the area so easy to compute!), is well defined even regardless of the condition number (so not an issue numerically?).
That is, plugin in 3.9.9 into 3.9.3 leads to a nicely well defined expression with no numerical issues. That seems like a plus for the "algebraic" form?
Notice that Nₐ,η and Nₐ,ξ are trivial for triangles leading to nice analytic expressions for the gradient. That is, gradient "times the jacobian (area)" seems to be very well defined.
Let me know if this does not make sense. I could expand Eq. 3.9.3 explicitly for you if that'd help.
There was a problem hiding this comment.
Reviewable status: 3 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy and @SeanCurtis-TRI)
geometry/proximity/volume_mesh.h, line 365 at r10 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
I'd recommend this instead:
For an arbitrary V, AB x AC may point into or out of the tetrahedron. In some cases it points in, in some cases it points out; that is a property of the definition of the tetrahedron. However, we are mathematically guaranteed that the gradient will always point inward because (AB x AC) / (AB x AC)⋅AV cancels out that variation.
For a given V, we can know whether AB x AC should point inward or outward. If the tetrahedron is degenerate (but not quite zero area), the numerical calculation may not match that classification. Rounding error may end up dominating the calculation producing a value for AB x AC pointing in the wrong direction or AV's direction may be essentially a random direction. In these degenerate cases, we aren't guaranteed inward pointing gradients and must treat them as zero.
And this is where we would be obliged to use the right epsilon value.
I added explanation in more details combining with ideas from your comments.
I have reservation on the phrase "that is a property of the definition of the tetrahedron" because it sounds to me like a positive-oriented tetrahedron always have AB x AC point inwards, and a negative-oriented tetrahedron always have AB x AC points outwards. However, it also depends on the local index i of V. On the same tetrahedron, there will be two i's that have AB x AC inwards and two i's that have AB x AC outwards.
geometry/proximity/volume_mesh.h, line 368 at r10 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: I really don't like this name. It's a scalar, a vector, and a shape all in one. And you opt not to use it in the next line (the triple product uses
ABxAC_M). I think this name is trying to do everything and ends up doing nothing.Ultimately, it seems like you're trying to make the code pedagogical. That someone could look at any single line and understand the algorithm and larger mathematical principles. At the very least, I'd bring this closer in line with your documentation above and eliminate "parallelogram" on this one and "parallelepiped" from the volume. It doesn't make me 100% happy, but it matches your documentation and I can be content with that.
I only keep area_vector and signed_volume for the names of the variables now.
geometry/proximity/volume_mesh.h, line 374 at r10 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW Based on my previous note about the possible negativity of this quantity, the
absprobably doesn't contribute anything.
Per f2f discussion, signed_volume can be a negative number with large absolute value, so we need abs.
geometry/proximity/test/mesh_intersection_test.cc, line 685 at r9 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW You should probably not use
T-- in DrakeTtends to have some very strong associations with scalar type and life would be simpler if you avoided it. Although it might be ok as it only appears in these couple of lines. I'm content in merely expressing the passing thought and defer to you on what you want to do with it.
Done. Changed T to F.
geometry/proximity/test/mesh_intersection_test.cc, line 706 at r9 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW Worth a reminder here in the test that
5π/8is our hard-coded threshold?
Done. Added that three lines above.
geometry/proximity/test/volume_mesh_test.cc, line 103 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: As for testing
Equal, you've only tested one way for them to be unequal. There are lots of expected conditions that would invalidate equality; you might consider either TODO'ing those or DO'ing those. :)
Done. Added more cases for testing.
There was a problem hiding this comment.
Reviewed 5 of 5 files at r11.
Reviewable status: 6 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy and @SeanCurtis-TRI)
geometry/proximity/volume_mesh.h, line 365 at r10 (raw file):
Previously, DamrongGuoy (Damrong Guoy) wrote…
I added explanation in more details combining with ideas from your comments.
I have reservation on the phrase "that is a property of the definition of the tetrahedron" because it sounds to me like a positive-oriented tetrahedron always have AB x AC point inwards, and a negative-oriented tetrahedron always have AB x AC points outwards. However, it also depends on the local index i of V. On the same tetrahedron, there will be two i's that have AB x AC inwards and two i's that have AB x AC outwards.
Essentially everything you said was exactly what I meant as "the definition of the tetrahedron". :) It's not a universal definition, it's just how we've defined it. Your solution is complete, clear, and elegant.
geometry/proximity/volume_mesh.h, line 378 at r11 (raw file):
// // If the tetrahedron is degenerate (tetrahedron with almost zero volume), // the calculation (AB x AC) / ((AB x AC)⋅AV) is not reliable any more due
BTW I'd maybe advocate the phrase "...is not numerically reliable...". But if you feel that's implicit in the word "calculation" I can accept that.
geometry/proximity/volume_mesh.h, line 392 at r11 (raw file):
constexpr double kEps3 = kEps * kEps * kEps; using std::abs; if (abs(signed_volume) <= kEps3) {
nit: We still need to account for the fact that a large degenerate triangle can still foil this test. As the volume gets larger, the scale of the rounding error increases. If the rounding error is scaled sufficiently large, it can pass this test but still have the wrong sign.
Per our f2f, you can indicate that we currently assume that our tets are generally much smaller than "unit size". However, if we get large degenerate triangles, this test won't be sufficient.
geometry/proximity/test/volume_mesh_test.cc, line 112 at r11 (raw file):
{ const auto mesh_translate = TestVolumeMesh<T>(RigidTransform<T>(Vector3<T>(1., 1., 1.)));
BTW If you wanted to strongly emphasize the fact that this test depends on bit exactitude in the values, you could do this use std::nextafter to change a single vertex position measure by a single bit.
geometry/proximity/test/volume_mesh_test.cc, line 137 at r11 (raw file):
tetrahedra.pop_back(); // 1 instead of 2 tetrahedra. VolumeMesh<T> another_mesh(std::move(tetrahedra), std::move(vertices)); EXPECT_FALSE(mesh->Equal(another_mesh));
nit: This single test doesn't test both different number of tets and vertices. It will only test different number of tets. The number of vertices are irrelevant to this test.
They should be broken into two different cases.
DamrongGuoy
force-pushed
the
field_gradient
branch
from
8eb59ad
to
c3b9ad7
Compare
There was a problem hiding this comment.
Reviewable status: 3 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @SeanCurtis-TRI)
a discussion (no related file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
Make sure there's a test on
mesh_intersectionthat the contact surface's gradients are reported in the world frame.
Done in (MeshIntersectionTest, ComputeContactSurfaceSoftRigidMoving). Now X_WS is no longer an identity transform, and we check the contact surface's gradients in World frame.
geometry/proximity/volume_mesh.h, line 365 at r10 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
Essentially everything you said was exactly what I meant as "the definition of the tetrahedron". :) It's not a universal definition, it's just how we've defined it. Your solution is complete, clear, and elegant.
Could you please unblock it now?
geometry/proximity/volume_mesh.h, line 378 at r11 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW I'd maybe advocate the phrase "...is not numerically reliable...". But if you feel that's implicit in the word "calculation" I can accept that.
Done. Added "numerically".
geometry/proximity/volume_mesh.h, line 392 at r11 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: We still need to account for the fact that a large degenerate triangle can still foil this test. As the volume gets larger, the scale of the rounding error increases. If the rounding error is scaled sufficiently large, it can pass this test but still have the wrong sign.
Per our f2f, you can indicate that we currently assume that our tets are generally much smaller than "unit size". However, if we get large degenerate triangles, this test won't be sufficient.
I put in a TODO and changed to kEps.
geometry/proximity/test/mesh_intersection_test.cc, line 408 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
So, the argument here is that this mesh field is only used in one test:
(MeshIntersectionTest, SampleVolumeFieldOnSurface)andSampleVolumeFieldOnSurface()doesn't do any math. The math is all done byClipTriangleByTetrahedron()andMeshFieldLinear::EvaluateCartesian(). And those functions are tested elsewhere. So, this doesn't need complex math, it just needs to exercise and poke.However, what this function does do is make sure that it's evaluating the field on the right tet and associating it with the right vertex. However, in this test you can't confirm that -- the nature of the intersection and pressure field means that for the three faces (and four vertices), the pressure is constant across all of them.
Done. I changed p4 to be 1000 times larger than p3, so the top tet has much lower pressure than the bottom tet. If SimpleVolumeFieldOnSurface picked the wrong tet, the test would fail.
geometry/proximity/test/volume_mesh_test.cc, line 137 at r11 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: This single test doesn't test both different number of tets and vertices. It will only test different number of tets. The number of vertices are irrelevant to this test.
They should be broken into two different cases.
Done. Two cases now.
There was a problem hiding this comment.
Reviewed 1 of 3 files at r12.
Reviewable status: 2 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @SeanCurtis-TRI)
geometry/proximity/test/mesh_intersection_test.cc, line 408 at r3 (raw file):
Previously, DamrongGuoy (Damrong Guoy) wrote…
Done. I changed p4 to be 1000 times larger than p3, so the top tet has much lower pressure than the bottom tet. If
SimpleVolumeFieldOnSurfacepicked the wrong tet, the test would fail.
That certainly distinguishes tets. :)
The other half is the fact that vertices 0, 1, and 2 have zero-values. That means the test below that makes use of this volume field ((MeshIntersectionTest, SampleVolumeFieldOnSurface)), cannot assess the unique contributions of those field values. We know the final answer must depend on all of them, generally. But in this case, weighted sums of zero look just like not doing addition. Or using the wrong weights, etc.
I know it makes articulating the test harder (it's not quite so nice to compute the expected answer by hand). But it means in this one test, all of the bits matter.
There was a problem hiding this comment.
Reviewed 2 of 3 files at r12.
Reviewable status: 8 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy)
geometry/proximity/test/mesh_intersection_test.cc, line 972 at r12 (raw file):
typename Mesh::VertexIndex* vertex) { for (typename Mesh::VertexIndex v(0); v < mesh_M.num_vertices(); ++v) { if ((p_MQ - mesh_M.vertex(v).r_MV()).norm() < 1e-14) {
BTW If you wanted to bump performance you could compare squared distance with 1e-28.
geometry/proximity/test/mesh_intersection_test.cc, line 998 at r12 (raw file):
// true if found. // @pre Assume Q is numerically well inside a tetrahedron or well outside the // mesh. It may incorrectly classify Q near vertices, edges, or faces of a
Is this still true with your allowing that any barycentric coordinate can be a little negative?
geometry/proximity/test/mesh_intersection_test.cc, line 1046 at r12 (raw file):
const auto X_WR = X_WS * X_SR; // Contact surface C is expressed in World frame. const auto contact_C_W = ComputeContactSurfaceFromSoftVolumeRigidSurface(
nit: This name is misleading: you have no frame C in this test and, by definition, contact surfaces have their mesh (and related quantities) define din the world frame (as indicated by the ContactSurface::mesh_W() method. So, contact_C_W is misleading because it's suggestive a non-existant frame C. And contact_W is likewise misleading because it suggests that you could have a contact_S = X_SW * contact_W (suggesting that the contact surface's mesh is measured and expressed in frame S which contradicts the API of the ContactSurface).
You might as well just call it contact_surface and be done with it.
This spills down below to other notation.
geometry/proximity/test/mesh_intersection_test.cc, line 1069 at r12 (raw file):
} // Check that the pressure gradient on the contact surface C equals
nit Here's the big question. Why are we computing the pressure gradient on the contact surface? I thought we weren't going to do that.
geometry/proximity/test/mesh_intersection_test.cc, line 1074 at r12 (raw file):
// TODO(DamrongGuoy): Remove this dynamic cast when we remove MeshField // and keep only MeshFieldLinear. auto pressure_C_W =
nit const auto* pressure_W is more suggestive of the facts.
Also, as noted above, the use of C in this name is somewhere between misleading and incorrect.
geometry/proximity/test/mesh_intersection_test.cc, line 1091 at r12 (raw file):
ASSERT_TRUE(FindElement(p_SQ, soft_S->mesh(), &tetrahedron, &b_Q)); // Pressure gradient vector in S expressed in S's frame. const Vector3d gradp_S_S = soft_S->EvaluateGradient(tetrahedron);
nit: This is another case where the double label is misleading. The _S_S suggests all the wrong things. You should limit the use of _S kinds of notation to actually mean frames. This might be: grad_soft_S (later transformed to grad_soft_W)? Or some such thing. As it is, you've implicitly equated p_S with soft without particular value.
geometry/proximity/test/mesh_intersection_test.cc, line 1098 at r12 (raw file):
// subtracting the normal component, i.e., ∇p_S - (∇p_S⋅n)n, where n is // the unit face normal of face 0 of C. const Vector3d nhat_C_W =
nit (And to complete the pedantic tour) nhat_C_W has similar issues.
geometry/proximity/test/mesh_intersection_test.cc, line 1125 at r12 (raw file):
// -Z // // To avoid "double counting" problem, we then translate R half a unit
nit You need an article. Either '...avoid the "double counting"...' or '...avoid a "double counting"...' (probably prefer the first).
DamrongGuoy
force-pushed
the
field_gradient
branch
from
c3b9ad7
to
cde3934
Compare
There was a problem hiding this comment.
Reviewable status: 1 unresolved discussion, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @SeanCurtis-TRI)
geometry/proximity/test/mesh_intersection_test.cc, line 408 at r3 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
That certainly distinguishes tets. :)
The other half is the fact that vertices 0, 1, and 2 have zero-values. That means the test below that makes use of this volume field (
(MeshIntersectionTest, SampleVolumeFieldOnSurface)), cannot assess the unique contributions of those field values. We know the final answer must depend on all of them, generally. But in this case, weighted sums of zero look just like not doing addition. Or using the wrong weights, etc.I know it makes articulating the test harder (it's not quite so nice to compute the expected answer by hand). But it means in this one test, all of the bits matter.
Done. Added this note in the test of SampleVolumeFieldOnSurface().
// Here we exploit the simplicity of TrivialVolumeMeshField<>() to check
// the field value. The test of field evaluation with more complex field
// values are in the unit test of VolumeMeshFieldLinear<>. Furthermore,
// testing that the right vertex is assigned the right field value is done
// in testing ComputeContactSurfaceFromSoftVolumeRigidSurface().geometry/proximity/test/mesh_intersection_test.cc, line 972 at r12 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW If you wanted to bump performance you could compare squared distance with 1e-28.
Ok.
geometry/proximity/test/mesh_intersection_test.cc, line 998 at r12 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
Is this still true with your allowing that any barycentric coordinate can be a little negative?
Ok. I forgot to change it when I flipped the sign of the tolerance.
geometry/proximity/test/mesh_intersection_test.cc, line 1046 at r12 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: This name is misleading: you have no frame
Cin this test and, by definition, contact surfaces have their mesh (and related quantities) define din the world frame (as indicated by theContactSurface::mesh_W()method. So,contact_C_Wis misleading because it's suggestive a non-existant frameC. Andcontact_Wis likewise misleading because it suggests that you could have acontact_S = X_SW * contact_W(suggesting that the contact surface's mesh is measured and expressed in frame S which contradicts the API of theContactSurface).You might as well just call it
contact_surfaceand be done with it.This spills down below to other notation.
Ok. Just contact now.
geometry/proximity/test/mesh_intersection_test.cc, line 1069 at r12 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit Here's the big question. Why are we computing the pressure gradient on the contact surface? I thought we weren't going to do that.
Ok. We'll fix it in the future PR.
geometry/proximity/test/mesh_intersection_test.cc, line 1074 at r12 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit
const auto* pressure_Wis more suggestive of the facts.Also, as noted above, the use of
Cin this name is somewhere between misleading and incorrect.
Done. const auto* c_pressure_W now.
geometry/proximity/test/mesh_intersection_test.cc, line 1091 at r12 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: This is another case where the double label is misleading. The
_S_Ssuggests all the wrong things. You should limit the use of_Skinds of notation to actually mean frames. This might be:grad_soft_S(later transformed tograd_soft_W)? Or some such thing. As it is, you've implicitly equatedp_Swithsoftwithout particular value.
Done. grad_s_pressure_S now.
geometry/proximity/test/mesh_intersection_test.cc, line 1098 at r12 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit (And to complete the pedantic tour)
nhat_C_Whas similar issues.
Done. nhat_W now.
geometry/proximity/test/mesh_intersection_test.cc, line 1125 at r12 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit You need an article. Either '...avoid the "double counting"...' or '...avoid a "double counting"...' (probably prefer the first).
Done.
There was a problem hiding this comment.
There's a note down below that I'll echo here since it's a fundamental enough issue. I thought we'd agreed that MeshFieldLinear::CalcGradientField() would become public, it wouldn't be automatically called in the constructor, and that instantiators of those would have to explicitly opt in.
The PR doesn't reflect that conclusion? Have I forgotten something? Or is there some other reason we're still doing the automatic generation?
Reviewed 1 of 1 files at r13.
Reviewable status: 2 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @DamrongGuoy)
geometry/proximity/test/mesh_intersection_test.cc, line 1069 at r12 (raw file):
Previously, DamrongGuoy (Damrong Guoy) wrote…
Ok. We'll fix it in the future PR.
Still confused. I thought our conversation was that the pressure gradient wouldn't be computed for a mesh automatically. One had to explicitly call it to compute it (and we simply wouldn't do that for ContactSurfaces. But I note it is sill being called automatically in the constructor.
Was your understanding of our resolution different?
geometry/proximity/test/mesh_intersection_test.cc, line 1027 at r13 (raw file):
// comprehensive way. GTEST_TEST(MeshIntersectionTest, ComputeContactSurfaceSoftRigidMoving) { // Soft octahedron volume S with pressure field.
nit: You are still using S to mean the mesh and not the frame of the mesh. (Similarly for R.) That's slightly unfortunate. The only relic that still relies on that is id_S (and, of course, id_R). If you continued the pattern, they'd become s_id and r_id, respectively and then S and R would unambiguously just be frames.
I'd suggest you go through all of the R and S instances and make sure that they are either properly referring to frames or change the documentation.
geometry/proximity/test/mesh_intersection_test.cc, line 1057 at r13 (raw file):
EXPECT_EQ(12, contact->mesh_W().num_faces()); // Point Q is at the center vertex of S. Check that C also has a vertex
nit: This is another confusing artifact due to the overload of S. In this case, you mean the mesh and not the frame.
DamrongGuoy
force-pushed
the
field_gradient
branch
from
cde3934
to
fee1538
Compare
There was a problem hiding this comment.
Done. Per f2f discussion, now we have a parameter in the constructor of MeshFieldLinear to control whether to calculate the gradient field. We also throw in EvaluateGradient() if the gradient was not calculated.
Reviewable status: LGTM missing from assignee SeanCurtis-TRI(platform), needs at least two assigned reviewers (waiting on @SeanCurtis-TRI)
geometry/proximity/test/mesh_intersection_test.cc, line 1069 at r12 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
Still confused. I thought our conversation was that the pressure gradient wouldn't be computed for a mesh automatically. One had to explicitly call it to compute it (and we simply wouldn't do that for
ContactSurfaces. But I note it is sill being called automatically in the constructor.Was your understanding of our resolution different?
Done. Removed it now.
geometry/proximity/test/mesh_intersection_test.cc, line 1027 at r13 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: You are still using
Sto mean the mesh and not the frame of the mesh. (Similarly forR.) That's slightly unfortunate. The only relic that still relies on that isid_S(and, of course,id_R). If you continued the pattern, they'd becomes_idandr_id, respectively and thenSandRwould unambiguously just be frames.I'd suggest you go through all of the
RandSinstances and make sure that they are either properly referring to frames or change the documentation.
Ok. s_id and r_id now. In the comment I still would like to give the name S to the soft geometry. S is both the name of geometry and also the name of its frame.
geometry/proximity/test/mesh_intersection_test.cc, line 1057 at r13 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: This is another confusing artifact due to the overload of
S. In this case, you mean the mesh and not the frame.
Done. soft mesh s_mesh_S now.
There was a problem hiding this comment.
You have a monster PR. You'll have to persuade platform reviewers to get on board or figure out a way to chop this in half.
Reviewed 3 of 3 files at r14.
Reviewable status: 2 unresolved discussions, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @DamrongGuoy)
geometry/proximity/mesh_field_linear.h, line 129 at r14 (raw file):
Vector3<T> EvaluateGradient( typename MeshType::ElementIndex e) const { if (e >= gradients_.size()) {
BTW This condition and message are not well aligned. I could properly compute the gradient but then pass in a bad element index and get the misleading error that gradients haven't been calculated. Probably:
if (gradients_.size() != mesh_W().num_elements) {
(or some such thing) would be a better test. Or, if you assume that you only compute a mesh field for a non-empty mesh, then
if (gradients_.size() == 0) {
would do it.
geometry/proximity/mesh_intersection.h, line 563 at r14 (raw file):
*surface_MN_M = std::make_unique<SurfaceMesh<T>>( std::move(surface_faces), std::move(surface_vertices_M)); bool calculate_gradient = false;
nit const (here and below).
geometry/proximity/test/mesh_intersection_test.cc, line 1057 at r13 (raw file):
Previously, DamrongGuoy (Damrong Guoy) wrote…
Done.
soft mesh s_mesh_Snow.
I'm marking myself as satisfied. But your solution is a bit overkill. soft mesh s_mesh_S is inherently redundant. Both soft mesh and s_mesh_S independently carry all the information. Both together deliver it twice. I just didn't like "of S" (as "S" is itself designed to be ambiguous).
geometry/proximity/test/mesh_intersection_test.cc, line 1058 at r14 (raw file):
// Point Q is at the center vertex of the soft mesh s_mesh_S. Check that // the contact surface C also has a vertex coincident with Q with the same
BTW Are we still calling the contact surface C?
|
Github shows 1256 lines. That's under the 1500 ceiling. |
DamrongGuoy
force-pushed
the
field_gradient
branch
from
5b0e4c6
to
6506a98
Compare
There was a problem hiding this comment.
Reviewable status: needs at least two assigned reviewers (waiting on @SeanCurtis-TRI)
geometry/proximity/mesh_field_linear.h, line 129 at r14 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW This condition and message are not well aligned. I could properly compute the gradient but then pass in a bad element index and get the misleading error that gradients haven't been calculated. Probably:
if (gradients_.size() != mesh_W().num_elements) {(or some such thing) would be a better test. Or, if you assume that you only compute a mesh field for a non-empty mesh, then
if (gradients_.size() == 0) {would do it.
Done. Changed to if (gradients_.size() == 0)
geometry/proximity/mesh_intersection.h, line 563 at r14 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit
const(here and below).
Done.
geometry/proximity/test/mesh_intersection_test.cc, line 1058 at r14 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
BTW Are we still calling the contact surface
C?
Yes. S for soft, R for rigid, C for contact.
There was a problem hiding this comment.
@sammy-tri for platform review please.
Reviewable status: needs at least two assigned reviewers (waiting on @SeanCurtis-TRI)
There was a problem hiding this comment.
+@sammy-tri for platform review please.
Reviewable status: LGTM missing from assignee sammy-tri(platform) (waiting on @sammy-tri and @SeanCurtis-TRI)
There was a problem hiding this comment.
Reviewed 2 of 2 files at r15.
Reviewable status: 1 unresolved discussion, LGTM missing from assignee sammy-tri(platform) (waiting on @DamrongGuoy and @sammy-tri)
geometry/proximity/mesh_field_linear.h, line 125 at r15 (raw file):
The gradient is a vector in R³ expressed in frame M. For surface meshes, it will particularly lie parallel to the plane of the corresponding triangle. @throw std::runtime_error if the gradient vector was not calculated.
nit: The "throwing" behavior should be documented and tested.
DamrongGuoy
force-pushed
the
field_gradient
branch
from
6506a98
to
28377e2
Compare
There was a problem hiding this comment.
Reviewable status: LGTM missing from assignee sammy-tri(platform) (waiting on @sammy-tri and @SeanCurtis-TRI)
geometry/proximity/mesh_field_linear.h, line 125 at r15 (raw file):
Previously, SeanCurtis-TRI (Sean Curtis) wrote…
nit: The "throwing" behavior should be documented and tested.
Done in mesh_field_linear_test.cc, GTEST_TEST(MeshFieldLinearTest, TestEvaluateGradientThrow).
There was a problem hiding this comment.
Reviewed 1 of 1 files at r16.
Reviewable status: LGTM missing from assignee sammy-tri(platform) (waiting on @sammy-tri)
There was a problem hiding this comment.
First pass complete, mainly just a doxygen question. Full disclosure, this was mostly a style/flow review, I did not attempt to fully understand the implementation.
Reviewed 1 of 24 files at r1, 2 of 7 files at r4, 2 of 17 files at r5, 5 of 12 files at r9, 1 of 5 files at r11, 2 of 3 files at r12, 1 of 3 files at r14, 2 of 2 files at r15, 1 of 1 files at r16.
Reviewable status: 1 unresolved discussion, LGTM missing from assignee sammy-tri(platform) (waiting on @DamrongGuoy)
geometry/proximity/penetration_doxygen.h, line 16 at r16 (raw file):
@{ @defgroup module_contact_surface Contact Surface
I don't see any content in the "Contact Surface" section when I click through the doxygen. Is this deliberate/known? If so, can we add a placeholder for now saying that the content is still under development to make clear the blank section isn't a bug?
DamrongGuoy
force-pushed
the
field_gradient
branch
from
28377e2
to
d8e9256
Compare
There was a problem hiding this comment.
Reviewable status: 1 unresolved discussion, LGTM missing from assignee sammy-tri(platform) (waiting on @sammy-tri and @SeanCurtis-TRI)
geometry/proximity/penetration_doxygen.h, line 16 at r16 (raw file):
Previously, sammy-tri (Sam Creasey) wrote…
I don't see any content in the "Contact Surface" section when I click through the doxygen. Is this deliberate/known? If so, can we add a placeholder for now saying that the content is still under development to make clear the blank section isn't a bug?
Done. Thank you. I added the new file contact_surface_doxygen.h as a placeholder.
Previously I had contact_surface_doxygen.h with very long content, but to shorten this PR, I decided to remove it and will work on it in another PR.
I do like the idea to have the file as a placeholder for now.
There was a problem hiding this comment.
Reviewed 1 of 1 files at r17.
Reviewable status: 1 unresolved discussion, LGTM missing from assignee sammy-tri(platform) (waiting on @sammy-tri)
There was a problem hiding this comment.
Reviewed 1 of 1 files at r17.
Reviewable status:complete! all discussions resolved, LGTM from assignees sammy-tri(platform),SeanCurtis-TRI(platform)
Resolve issue #12441 "Hydroelastic contact surface broken for thin rigid object -- needs to use normals".
The size of this PR is much larger than originally expected. To keep it manageable, another PR will have the code for a new scene_graph example that shows the soft ball touching the thin plate and bowl. This PR only has the code and the doxygen doc to support that scene_graph example.
This change is