We introduce the concept of shape optimization of singular points of a weak solution in boundary value problems of partial differential equations (BVP) including fracture mechanics, shape optimization of boundary and interface, etc. Roughly speaking, singularity is a gap between weak and strong solutions. We also introduce the integral formula named the GJ-integral defined with an arbitrary subdomain as an extension of the J-integral in fracture mechanics. The GJ-integral is defined locally, is made from partial differential equations in BVP only, takes zero if the solution is regular, and expresses energy shape sensitivity. A generalization of the Hadamard variational formula is proposed and various cost functions are rewritten using the GJ-integral. In the last section there are numerical calculations using the GJ-integral with the finite element method.