Strain-gradient-driven electric polarization in soft dielectric polymers is especially pronounced at small scales. In such materials, size effects emerge not only from strain gradients but also from surface contributions, driven by their high surface-to-volume ratios and the increased role of surface tension. Capturing these effects requires theoretical and numerical frameworks that extend beyond classical elasticity. In this work, we present a unified isogeometric framework for modeling size-dependent mechanics in soft materials using two approaches: strain-gradient hyperelasticity and surface elasticity theory. The isogeometric formulation provides the necessary $C^1$ continuity to support higher-order spatial derivatives in the discretization essential for modeling strain-gradient-induced polarization, as well as accurate geometric and kinematic representation of material surfaces. This models a coupled boundary and bulk formulation under large strain assumptions, where a zero-thickness shell is perfectly bonded to a bulk material, where this curvature-resisting surface further necessitates a $C^1$-continuous discretization of the surface. We demonstrate this framework over examples of highly deformable soft polymers and biomimetic systems undergoing large deformations and instabilities such as buckling and wrinkling. In soft polymers, we compute strain-gradient drivent electric potential differences generated during deformation—linking mechanics to functional electromechanical response. Numerical simulations demonstrate the emergence of strong boundary-layer stresses and size-dependent instability patterns. These results underscore the importance of higher-order and surface effects in soft material systems and highlight the advantages of isogeometric analysis in resolving the size-dependent mechanics and electrostatics of soft solids.