Structures of faults, fractures and geomaterial interfaces are complex, making their behaviour difficult to predict. The ability to simulate geomechanical, geochemical and hydrogeological behaviours with these structures is required for better engineering process ahead of design. A pivotal step of geological simulation using the commonly applied finite element method is the mesh generation which is the subdivision of spatial objects into small elements. When complicated constraints such as faults, fractures and interfaces are present existing mesh generation methods are still far from feasible solutions for geological modelling. This thesis addresses this key issue and demonstrates how to generate high-quality mesh for representation of geo-objects with heterogeneous structures. This thesis focuses on proposing a comprehensive solution for meshing and remeshing for geological modelling. The primary contributions are summarized as follows:1. A boundary focused quadrilateral mesh generation method has been developed to produce highquality mesh models for 2D multi-material geological data. Together with geodesic isolines, a valence clear-up pattern "Pisces" developed in this thesis is employed to improve the mesh quality.Hence the proposed method has more advantages in generating high-quality elements compared with algorithms developed by previous researchers (Park et al. 2007, Lee et al. 2003 In the application to meshing images of coal plugs and fractured rocks, layers of well-aligned elements parallel to the material interfaces are generated and thus the average element quality of the generated mesh is close to that of a regular quadrilateral;2. An effective 3D meshing method has been proposed and implemented towards generating mesh representations for fractured rocks. A high-resolution image can capture fracture structures in detail but it leads to a huge data set. To reduce the data size, surface meshes (served as constraints within volume meshes) rather than traditional volume meshes are utilized to represent fractures in this research. A simplified Voronoi diagram based on a proposed pseudo-surface assumption is developed to extract fractures, identify fracture junctions and indirectly generate high-quality meshes. Subsequent numerical experiments showed that the generated mesh representing a fractured rock has a high shape similarity 64.57% and its data size is merely 1/6000 as much as that of the image data. II 3. A strategy has been carried out for meshing 3D geological reservoirs with arbitrary stratigraphic surface constraints. To achieve this, a geodesic-based procedure is designed to provide a robust implementation for handling geometrical heterogeneity on stratigraphic surfaces (i.e. stratum interfaces). The procedure also contributes to a well-designed advancing front technique which achieves the generation of high-quality elements and adaptive meshes. The effectiveness of the proposed strategy is demonstrated by meshing a reservoir geological model of Lawn Hill in Queensland, Australia. 4. A 3D reme...