The increasing geometric complexities in modern applications pose significant challenges to numerical analysis and computational science when using traditional techniques with either structured or unstructured meshes. We present our recent work towards a unified theoretical and algorithmic framework for accurate and stable numerical discretizations and efficient solution techniques for partial differential equations (PDEs) over complex geometries. This talk is composed of two parts. The first part will describe a unified and versatile framework of high-order accurate, stable, and efficient numerical methods based on local weighted-least squares (WLS) approximations and a global weighted-residual formulation for discretizing elliptic and parabolic PDEs over complex and curved geometries. This framework generalizes the finite difference methods to unstructured meshes in a coherent way, and generalizes the finite element methods to have relaxed requirements on mesh quality and to deliver higher-order convergence on curved geometries. The second part describes a hybrid geometric+algebraic multigrid method (or HyGA) for weighted-residual methods with hierarchical basis functions. HyGA combines the rigor, accuracy and runtime-and-memory efficiency of geometric multigrid with the robustness and flexibility of algebraic multigrid, and at the same time is relatively easy to implement. We also discuss the intricate interactions between the WLS-based PDE discretizations with multigrid solvers, including the structure of the linear systems and the solutions of rank-deficient linear systems.