$\newcommand{\Dev}{\operatorname{Dev}}$ $\newcommand{\dev}{\operatorname{dev}}$ $\newcommand{\grad}{\operatorname{grad}}$ $\newcommand{\Grad}{\operatorname{Grad}}$ $\newcommand{\divg}{\operatorname{div}}$ $\newcommand{\Ei}{\operatorname{Ei}}$ $\newcommand{\tr}{\operatorname{tr}}$ $\newcommand{\Divg}{\operatorname{Div}}$ $\newcommand{\cay}{\operatorname{cay}}$ $\newcommand{\rot}{\operatorname{rot}}$

The determination of the penalty factor $\varepsilon_{N}$ can be a difficult task, since a good value may depend on both material parameters and geometrical factors. In FEBio the value of this penalty factor can be determined automatically. In this case FEBio will calculate a penalty factor for each facet using the following formula. $$\begin{equation} \varepsilon_{i}=\frac{f_{SI}E_{i}A_{i}}{V_{i}}\,.\label{eq595} \end{equation}$$ Here, $E_{i}$ is the effective Young's modulus along the facet normal, $A_{i}$ the surface area of the facet, $V_{i}$ the volume of the element to which this facet belongs and $f_{SI}$ a user defined scale factor. The parameter $E_{i}$ is evaluated from the elasticity tensor $\boldsymbol{\mathcal{C}}$ and the facet unit normal $\mathbf{n}$ according to $$\begin{equation} \frac{1}{E_{i}}=\left(\mathbf{n}\otimes\mathbf{n}\right):\boldsymbol{\mathcal{C}}^{-1}:\left(\mathbf{n}\otimes\mathbf{n}\right)\,,\label{eq596} \end{equation}$$ where $\boldsymbol{\mathcal{C}}^{-1}$ is the compliance tensor.