Possibly inconsistent results for convergent formulation with respect to dhat?

[Andlon]

Possibly inconsistent results for convergent formulation with respect to dhat?

[Andlon]

Ah, that's an interesting interpretation! Here's what I've been thinking about the value of $\kappa$ (when not being automatically determined). It should:

• be high enough so that, at force equilibrium, for all distances, $d^2 / \hat{d}^2 > 0.05$ (or something like that, I can't remember what a reasonable value is), i.e. we want to keep the distance in the "well-conditioned" part of the log barrier (i.e. in the part of the domain where it does not change too rapidly).
• not be so high as to cause excessive stiffness and conditioning issues in the other direction

Since we look at the ratio between $d / \hat{d}$ the $\kappa$ value should be relatively independent of choice of $\hat{d}$ (and I've observed that in some admittedly too simple experiments). But I haven't tried to vary elastic parameters yet. I think I lack the right physical intuition to understand why material stiffness impacts the "right" value of $\kappa$. Intuitively I would think that - at least for resting contacts - the stiffness does not matter so much for the relative location of the force equilibrium. I have to think about this!

Feel free to close this unless you want to track the change you said you wanted to make. I guess long term you probably also want to document what the convergent formulation actually does, but I suspect this will eventually appear in due time anyway :-)

[zfergus]

The intention of this change is to be able to express the barrier potential in the same physical units as the elasticity potential (i.e., Pa $\cdot$ m) (more on this at the bottom). With non-squared distances, we do this by writing the barrier function as

$$ b_\text{new}(d, \hat{d}) = -\kappa\hat{d}\left(\frac{d}{\hat{d}} - 1\right)^2 \log\left(\frac{d}{\hat{d}}\right). $$

The relationship between $b_\text{old}(d, \hat{d})$ and $b_\text{new}(d, \hat{d})$ is

$$ b_\text{new}(d, \hat{d}) = \frac{1}{\hat{d}} b_\text{old}(d, \hat{d}). $$

In the toolkit, we only implement $b_\text{old}(d, \hat{d})$ and store this division in the weight for simplicity of implementation.

However, we are using distances squared so we actually want to compute

$$ b_\text{new}(d^2, \hat{d}^2) = -\kappa\hat{d}\left(\frac{d^2}{\hat{d}^2} - 1\right)^2 \log\left(\frac{d^2}{\hat{d}^2}\right). $$

Again the relationship between $b_\text{old}(d^2, \hat{d}^2)$ and $b_\text{new}(d^2, \hat{d}^2)$ is

$$ b_\text{new}(d, \hat{d}) = \frac{1}{\hat{d}^3} b_\text{old}(d, \hat{d}). $$

Again, we only implement $b_\text{old}(d, \hat{d})$, so we want to store this division in the weight for simplicity of implementation. However, I am dividing by $\hat{d}^2$ when I should be dividing by $\hat{d}^3$ to get the correct units. I will make this change.

This will still not produce the same behavior as you expect (because we multiply the normalized barrier by $\hat{d}$), but this is the correct implementation of the algorithm we want (i.e., physically-based units for the barrier term). Alternatively, if you wanted to use the normalized barrier, you could assume the units of $\kappa$ are Pa $\cdot$ m. What I like about the intended implementation is that $\kappa$ has units of Pa and we can express $\kappa$ as a ratio of the materials Young's modulus ($0.1E$ works well in many examples).

The intention is to treat the barriers as a thin elastic region around the mesh, and having consistent units makes it easier to pick the stiffness for this "material".