A formulation based on the approximation of the stress field is used to compute directly the stress intensity factors in crack problems. The boundary displacements are independently approximated. In each finite element, the assumed stresses may model multipoint singularities of variable order. The differential equilibrium equations are locally satisfied as solutions of the governing differential system are used to build the stress approximation basis. The approximation on the boundary displacements is constrained to satisfy locally the kinematic boundary conditions. The remaining fundamental conditions, namely the differential compatibility equations, the constitutive relations and the static boundary conditions are enforced through weighted residual statements. The approximation criteria are so chosen as to ensure that the finite element model is described by a sparse, adaptive and symmetric governing system described by structural matrices with boundary integral expressions. Numerical applications are presented to show that accurate solutions can be obtained using structural discretizations based on coarse meshes of few but highly rich elements, each of which may have different geometries and alternative approximation laws.