In many problems of geophysical interest, one has to deal with data that exhibit complex fault structures. This occurs, for instance, when describing the topography of seafloor surfaces, mountain ranges, volcanoes, islands, or the shape of geological entities, as well as when dealing with reservoir characterization and modelling. In all these circumstances, due to the presence of large and rapid variations in the data, attempting a fitting using conventional approximation methods necessarily leads to instability 68 Numer Algor (2008) 48:67-92 phenomena or undesirable oscillations which can locally and even globally hinder the approximation. As will be shown in this paper, the right approach to get a good approximant consists, in effect, in applying first a segmentation process to precisely define the locations of large variations and faults, and exploiting then a discrete approximation technique. To perform the segmentation step, we propose a quasi-automatic algorithm that uses a level set method to obtain from the given (gridded or scattered) Lagrange data several patches delimited by large gradients (or faults). Then, with the knowledge of the location of the discontinuities of the surface, we generate a triangular mesh (which takes into account the identified set of discontinuities) on which a D m -spline approximant is constructed. To show the efficiency of this technique, we will present the results obtained by its application to synthetic datasets as well as real gridded datasets in Oceanography and Geosciences.