In　this　study,　the　two-dimensional　physical　domain　containing　cracks　is　divided　into　several　non　overlapping　parts:　rectangular　crack-tip　regions　around　crack　tips　and　the　outer　region　without　any　crack　tip.　In　each　crack-tip　region　the　displacement　is　approximated　with　Williams＇　series;　while　in　the　outer　region　it　is　approximated　with　numerical　manifold　interpolation.　In　order　to　balance　accuracy　and　efficiency　in　solution,　a　transitional　zone　encompassing　each　crack-tip　region　is　locally　refined　with　a　structured　mesh.　To　avoid　singular　integration　over　a　crack-tip　region,　the　potential　energy　over　every　crack-tip　region　is　transformed　into　the　boundary　integration.　Three　different　methods　to　enforce　compatibility　on　interfaces　are　compared,　concluding　the　Lagrange　multiplier　method　is　superior　over　the　other　two.　(C)　2016　Elsevier　Ltd.　All　rights　reserved.