Physics-Based Simulation

Solid-Obstacle Contact

Recall that we used a conservative force model to approximate the contact traction ${\mathbf{T}}_C$, allowing it to be directly evaluated given the current configuration of the solids. This results in a contact potential:

$$P_C={\int }_{{\Gamma }_C}\frac{1}{2}b\left(\underset{{\mathbf{X}}_2\in {\Gamma }_C-\mathcal{N}(\mathbf{X})}{\min }\Vert \mathbf{x}(\mathbf{X},t)-\mathbf{x}({\mathbf{X}}_2,t)\Vert ,\hat{d}\right)ds(\mathbf{X}),$$

where $b()$ is the barrier energy density function, and $\mathcal{N}(\mathbf{X})$ is an infinitesimal region around $\mathbf{X}$ where contact is ignored for theoretical soundness.

For normal contact between simulated solids and collision obstacles (ignoring self-contact for now), $P_C$ can be written in a much simpler form $$\begin{array}{rl}{P}_C& ={\int }_{{\Gamma }_S}\frac{1}{2}b(\underset{{\mathbf{X}}_2\in {\Gamma }_O}{\min }\Vert \mathbf{x}(\mathbf{X},t)-\mathbf{x}({\mathbf{X}}_2,t)\Vert ,\hat{d})ds(\mathbf{X})\\ & +{\int }_{{\Gamma }_O}\frac{1}{2}b(\underset{{\mathbf{X}}_2\in {\Gamma }_S}{\min }\Vert \mathbf{x}(\mathbf{X},t)-\mathbf{x}({\mathbf{X}}_2,t)\Vert ,\hat{d})ds(\mathbf{X})\\ & ={\int }_{{\Gamma }_S}b(\underset{{\mathbf{X}}_2\in {\Gamma }_O}{\min }\Vert \mathbf{x}(\mathbf{X},t)-\mathbf{x}({\mathbf{X}}_2,t)\Vert ,\hat{d})ds(\mathbf{X})\\ & ={\int }_{{\Gamma }_S}b(d^{PO}(\mathbf{x}(\mathbf{X},t),O),\hat{d})ds(\mathbf{X}).\end{array}$$ Here ${\Gamma }_S$ and ${\Gamma }_O$ are the boundaries of the simulated solids and obstacles respectively, $d^{PO}(\mathbf{x}(\mathbf{X},t),O)={\min }_{{\mathbf{X}}_2\in {\Gamma }_O}\Vert \mathbf{x}(\mathbf{X},t)-\mathbf{x}({\mathbf{X}}_2,t)\Vert$ is the point-obstacle distance, and the simplification from two terms to one single term is due to symmetry in the continuous setting. With triangle discretization, $${\int }_{{\Gamma }_S}b(d^{PO}(\mathbf{x}(\mathbf{X},t),O),\hat{d})ds(\mathbf{X})\approx \sum _{e\in \mathcal{T}}{\int }_{\partial {\Omega }_e^0\cap {\Gamma }_S}b(d^{PO}(\mathbf{x}(\mathbf{X},t),O),\hat{d})ds(\mathbf{X}).\text{\hspace{1em}(20.2.1)}$$ Similar to the derivation for Neumann boundaries, for any boundary node $\hat{a}$, with 2 incident triangles, let us look at one of the integral. Without loss of generality, we can assume $N_{\hat{a}}=\beta$ (${\mathbf{X}}_{\hat{a}}$ corresponds to ${\mathbf{X}}_2$ in triangle $e$), and that $X_3$ is the other node of $e$ on the boundary edge. Then, switching the integration variables to $\beta$ gives us $$\begin{array}{rl}& {\int }_{\partial {\Omega }_e^0\cap {\Gamma }_S}b(d^{PO}(\mathbf{x}(\mathbf{X},t),O),\hat{d})ds(\mathbf{X})\\ =& {\int }_0^{1}b(d^{PO}(\mathbf{x}(\beta {\mathbf{X}}_2+(1-\beta ){\mathbf{X}}_3,t),O),\hat{d})\left|\frac{\partial s}{\partial \beta }\right|d\beta .\end{array}\text{\hspace{1em}(20.2.2)}$$ Since $b()$ and $d^{PO}()$ are both highly nonlinear functions, we could not obtain a closed-form expression for Equation (20.2.2). If we take the two end points ${\mathbf{X}}_2$ and ${\mathbf{X}}_3$ as quadrature points both with weights $\frac{1}{2}$, we can approximate the integral as $$\begin{array}{rl}& {\int }_0^{1}b(d^{PO}(\mathbf{x}(\beta {\mathbf{X}}_2+(1-\beta ){\mathbf{X}}_3,t),O),\hat{d})\left|\frac{\partial s}{\partial \beta }\right|d\beta \\ \approx & \frac{1}{2}b(d^{PO}(\mathbf{x}({\mathbf{X}}_2,t),O),\hat{d})\left|\frac{\partial s}{\partial \beta }\right|+\frac{1}{2}b(d^{PO}(\mathbf{x}({\mathbf{X}}_3,t),O),\hat{d})\left|\frac{\partial s}{\partial \beta }\right|.\end{array}\text{\hspace{1em}(20.2.3)}$$ Then, the whole boundary integral can be approximated as $${\int }_{{\Gamma }_S}b(d^{PO}(\mathbf{x}(\mathbf{X},t),O),\hat{d})ds(\mathbf{X})\approx \sum _{\hat{a}}\frac{\Vert {\mathbf{X}}_{\hat{a}}-{\mathbf{X}}_{\hat{a}-1}\Vert +\Vert {\mathbf{X}}_{\hat{a}}-{\mathbf{X}}_{\hat{a}+1}\Vert }{2}b(d^{PO}({\mathbf{x}}_{\hat{a}},O),\hat{d}),$$ assuming that ${\mathbf{X}}_{\hat{a}-1}$ and ${\mathbf{X}}_{\hat{a}+1}$ are the two neighbors of ${\mathbf{X}}_{\hat{a}}$ on the boundary. This is now exactly what has been implemented in Filter Line Search.

Remark 20.2.1 (Quadrature Choice for Line Segment). Selecting the two end points ($\beta =0,1$) as quadrature points for a line segment integral (Equation (20.2.3)) is not a common design choice. Typically, Gaussian quadrature would use $\beta =\frac{3\pm \sqrt{3}}{6}$. The advantage of choosing $\beta =0,1$ is that it results in fewer quadrature points globally, thus reducing computational costs, as neighboring edges share end points.

To see how $P_C$ connects to the boundary integral (Equation (20.1)) in the discrete weak form, let us take the derivative of the discretized contact potential (Equation (20.2.1)) with respect to ${\mathbf{x}}_{\hat{a}}$: $$\begin{array}{rl}& -\frac{\partial ({\sum }_{e\in \mathcal{T}}{\int }_{\partial {\Omega }_e^0\cap {\Gamma }_S}b(d^{PO}(\mathbf{x}(\mathbf{X},t),O),\hat{d})ds(\mathbf{X}))}{\partial {\mathbf{x}}_{\hat{a}}}\\ =& \sum _{e\in \mathcal{T}}{\int }_{\partial {\Omega }_e^0\cap {\Gamma }_S}-\frac{\partial b(d^{PO}(\mathbf{x}(\mathbf{X},t),O),\hat{d})}{\partial \mathbf{x}}\frac{\partial \mathbf{x}}{\partial {\mathbf{x}}_{\hat{a}}}ds(\mathbf{X})\\ =& \sum _{e\in \mathcal{T}}{\int }_{\partial {\Omega }_e^0\cap {\Gamma }_S}-\frac{\partial b(d^{PO}(\mathbf{x}(\mathbf{X},t),O),\hat{d})}{\partial \mathbf{x}}{N}_{\hat{a}}(\mathbf{X})ds(\mathbf{X}).\end{array}$$ Then we also verified that ${\mathbf{T}}_C(\mathbf{X},t)=-\frac{\partial b(d^{PO}(\mathbf{x}(\mathbf{X},t),O),\hat{d})}{\partial \mathbf{x}}$ here.