In recent years, Very Large Floating Structures (VLFSs) such as floating bridges and Mega-floats have been planned for efficient utilization of ocean space. Accurate prediction of elastic behavior of VLFSs due to wave action is indispensable for their design owing to their flexibility. The present study has been carried out under the Mega-float Project as one for developing an accurate general-purpose three-dimensional computer program for hyroelastic responses of VFLSs due to wave action in the water area with complicated geometry close to the real situation. This paper presents a numerical method for predicting the hydroelasic behavior of a VLFS in the complicated water area sheltered partially by breakwaters and land. A standard modal approach is adapted for the fluid-structure interactions. For the free-surface flow, an improved hybrid finite/infinite element formulation has been developed to reduce the computational size and cost for practical applications substantially. The flow field is divided by a fictitious boundary into an outer infinite domain and an inner irregular one, which are again divided into subdomains appropriately. The velocity potentials for the inner and outer subdomains are approximated by the corresponding vertical orthogonal eigenfunction expansions, and planar finite elements and hybrid infinite elements, respectively, and are matched each other on the corresponding interfaces in an averaged sense. To very large rectangular plates, analytical Rayleigh-Ritz methods are applied. Numerical results for a VLFS in the open sea have shown a satisfactory agreement with the existing numerical and measured ones. The response of a VLFS in a shelterd area and wave diffraction/refraction have also been examined. Conclusions are that the present 3D hybrid element code is a powerful working tool for practical VLFS investigations. Advanced simulations for 3D VLFS models will be reported in the subsequent papers.