Three-dimensional structural modeling is gaining importance for a broad range of quantitative geoscientific applications. However, existing approaches are still limited by the type of structural data they are able to use and by their lack of structural meaning. Most techniques heavily rely on spatial data for modeling folded layers, but are unable to completely use cleavage and lineation information for constraining the shape of modeled folds. This lack of structural control is generally compensated by expert knowledge introduced in the form of additional interpretive data such as cross-sections and maps. With this approach, folds are explicitly designed by the user instead of being derived from data. This makes the resulting structures subjective and deterministic. This paper introduces a numerical framework for modeling folds and associated foliations from typical field data. In this framework, a parametric description of fold geometry is incorporated into the interpolation algorithm. This way the folded geometry is implicitly derived from observed data, while being controlled through structural parameters such as fold wavelength, amplitude and tightness. A fold coordinate system is used to support the numerical description of fold geometry and to modify the behavior of classical structural interpolators. This fold frame is constructed from fold-related structural elements such as axial foliations, intersection lineations, and vergence. Poly-deformed terranes are progressively modeled by successively modeling each folding event going backward through time. The proposed framework introduces a new modeling paradigm, which enables the building of three-dimensional geological models of complex poly-deformed terranes. It follows a process based on the structural geologist approach and is able to produce geomodels that honor both structural data and geological knowledge.