Bivariate cubic L1 splines can provide shape-preserving surfaces for various applications. Using the reduced Hsieh-Clough-Tocher (rHCT) elements on the triangulated irregular networks (TINs), we model a bivariate cubic L1 spline as the solution to a nonsmooth convex programming problem. This problem is a generalized geometric programming (GGP) problem, whose dual problem is to optimize a linear objective function over convex cubic constraints. Using a linear programming transformation, a dual optimal solution can be converted to a desired primal solution. For computational efficiency, we further develop a compressed primal-dual interior-point method to directly calculate an approximated primal optimal solution. This compressed primal-dual algorithm can handle terrain data over hundreds-by-hundreds grids using a personal computer. However, for real-life applications, terrain data are given in thousands-by-thousands grids. To meet the computational challenge, we establish a "non-iterative" domain decomposition principle to reduce the computational requirements. We have also conducted computational experiments to show that the proposed domain decomposition principle can handle large size data for real terrain applications.