Since　the　shape　functions　derived　from　the　partition　of　unity-based　meshless　method,　such　as　the　element-free　Galerkin　method,　are　free　of　the　Kronecker　delta　property,　there　are　great　troubles　in　the　exact　imposition　of　the　essential　boundary　condition　and　boundary　continuity　of　materials.　Nevertheless,　if　adopting　the　penalty　method　or　the　Lagrange　multiplier　method,　problems,　like　the　selection　of　proper　penalty　factor　and　the　satisfaction　of　the　inf-sup　condition,　will　occur.　This　study　utilizes　the　property　of　partition　of　unity　that　once　the　local　solutions　satisfy　some　condition,　the　global　solution　will　automatically　satisfy　the　same　condition.　By　constructing　local　approximations　in　physical　patches　of　different　types　according　to　the　boundary　condition,　a　new　moving　least　square　interpolation-based　numerical　manifold　method(MLS-NMM)　is　developed.　Through　the　solution　of　some　typical　seepage　flow　problems,　it　is　demonstrated　that　the　proposed　procedure　is　capable　to　deal　with　the　problems　of　the　singular　angular　point　precisely　and　may　provide　an　alternative　solution　for　the　seepage　analysis　in　engineering.　©　2017,　Editorial　Office　of　Chinese　Journal　of　Geotechnical　Engineering.　All　right　reserved.