This paper presents a novel, automatic, level-set-based approach to model three-dimensional boundary problems and generate a new conformal body-fitted mesh. The proposed methodology is applied to long-term in-flight ice accretion, characterized by the formation of extremely irregular ice shapes. The level-set function is defined across the computational volume by comparing the position of the generic volume node with respect to the Lagrangian deformed geometry. This method avoids mesh entanglements and grid intersections typical of mesh deformation techniques, making it suitable for generating a body-fitted discretization of arbitrarily complex geometries as in-flight ice shapes. Effective and tunable smoothing techniques are proposed for dealing with complex ice shapes. Multistep numerical simulations over a horizontal stabilizer with deflected elevator and over a NACA0012 swept wing, both in rime and glaze conditions, are presented. The latter two are also compared with the experimentally measured ice shapes from the first AIAA Ice Prediction Workshop.