<p>Natural and synthetic polycrystalline materials exhibit complex mechanical behaviors governed by their microstructural features and interfacial properties. This paper introduces a robust computational framework for generating and analyzing polycrystalline structures with cohesive zone modeling capabilities in complex non-convex domains with arbitrary topology. The proposed framework efficiently handles non-smooth, disconnected domains with internal holes or gaps by computing the clipped Voronoi diagram, enabling the generation of complex computational models with realistic grain morphologies. The framework also provides several efficient algorithms to optimize and adjust the polygonal grains, making them suitable for finite element simulations and seamless insertion of both zero-thickness and finite-thickness cohesive zones along grain boundaries. We demonstrate the framework’s versatility through challenging applications, such as modeling calcium plaques in human arteries, where accurate representation of grain-grain and grain-matrix interfaces is critical. To promote reproducibility and facilitate broader adoption, we provide an open-source Python implementation that can be easily integrated with the scripting interfaces of commercial finite element solvers, offering researchers a powerful tool for studying fracture mechanisms in complex biological and synthetic composites.</p>