## 1 Introduction  
Polycrystalline materials are composed of several grains with sizes that range from several nanometers to millimeters and typically have random crystallographic orientations [1]. Most inorganic materials, including all common metals, most ceramics, rocks, and ice, are polycrystalline in nature and are composed of a large number of crystallites joined together by grain boundaries [2]. The mechanical behavior of polycrystalline materials is highly dependent on their microstructure, specifically the size and orientation of grains, as well as the area of grain boundaries. For instance, the Hall-Petch relationship states that the strength of polycrystals increases as the average grain size decreases [3]. Additionally, the manufacturing process plays a crucial role in determining the microstructural features of polycrystalline materials. This includes factors such as the cooling rate and the extent of deformation, which can greatly impact macroscopic properties of the material such as strength and toughness. Various numerical techniques have been developed to study the micromechanical response of polycrystals, including mean-field approaches (e.g., Taylor and self-consistent schemes) [4, 5, 6], Crystal Plasticity Finite Element (CPFE) [7], and [PERSON]'s function Fast Fourier Transform (FFT) models [8, 9]. While CPFE and FFT models are more computationally demanding than mean-field approaches, they can accurately capture the evolution of individual grain shapes, the interaction between neighboring grains, intra-granular strain evolution, and the heterogeneous onset of damage [10]. CPFE models provide the advantage of incorporating complex microstructural features and heterogeneities, such as micro-cracks, inclusions, and precipitates, as well as enabling simulations with complex boundary conditions resulting from inter/intra-granular interactions [11; 12]. However, creating image-based finite element (FE) models can be challenging. It is crucial to accurately represent grain morphologies, including grain boundaries and corners, when creating a high-fidelity FE model for CPFE simulations. Low resolution or the presence of noise in pixelated/voxelated imaging data obtained from sources such as electron backscatter diffraction (EBSD) can greatly hinder the achievement of this goal. Virtually reconstructed microstructures using methods such as Voronoi tessellation [13; 14] or DREAM.3D [15] can lack accuracy or lead to similar challenges when generating high-fidelity FE models. Reconstruction algorithms relying on the Voronoi tessellation method approximate the grain boundaries using straight planar surfaces, which is not realistic for many complex polycrystalline microstructures. Additionally, DREAM.3D synthesizes grain models with jagged boundaries (similar to EBSD data), and does not capture their edge/corner geometry. To bypass these challenges, CPFE simulations typically use pixelated/voxelated meshes [13; 14]. However, these meshes cannot accurately represent stress concentrations at the grain boundaries, which can negatively impact the accuracy of simulations aimed at modeling damage/debonding along grain boundaries. Several meshing algorithms are implemented to generate high-quality conforming FE meshes, such as Delaunay triangulation [16; 17], Quadtree/Octree-based approaches [18; 19], and advancing front methods [18; 20; 21]. The Delaunay triangulation method solves a nonlinear system of equations to optimize the mesh node locations, with the objective of ensuring that no mesh node is located within the circumcircle of a triangular element in 2D or the circumference of a tetrahedral element in 3D. In the advancing front method, new triangular/tetrahedral elements are created on an active front that emanates from specific regions in the domain, usually the boundaries, and gradually advance towards the interior to discretize the whole domain.