  Computing equilibrium shapes of crystals (ESC) is a challenging problem in
materials science that involves minimizing an orientation-dependent (i.e.,
anisotropic) surface energy functional subject to a prescribed mass constraint.
The highly nonlinear and singular anisotropic terms in the problem make it very
challenging from both the analytical and numerical aspects. Especially, when
the strength of anisotropy is very strong (i.e., strongly anisotropic cases),
the ESC will form some singular, sharp corners even if the surface energy
function is smooth. Traditional numerical approaches, such as the $H^{-1}$
gradient flow, are unable to produce true sharp corners due to the necessary
addition of a high-order regularization term that penalizes sharp corners and
rounds them off. In this paper, we propose a new numerical method based on the
Davis-Yin splitting (DYS) optimization algorithm to predict the ESC instead of
using gradient flow approaches. We discretize the infinite-dimensional
phase-field energy functional in the absence of regularization terms and
transform it into a finite-dimensional constraint minimization problem. The
resulting optimization problem is solved using the DYS method which
automatically guarantees the mass-conservation and bound-preserving properties.
We also prove the global convergence of the proposed algorithm. These desired
properties are numerically observed. In particular, the proposed method can
produce real sharp corners with satisfactory accuracy. Finally, we present
numerous numerical results to demonstrate that the ESC can be well simulated
under different types of anisotropic surface energies, which also confirms the
effectiveness and efficiency of the proposed method.
