The result
[tex]S = \frac{c}{3} \log \frac{\ell}{a}[/tex]
only applies to (1+1)-d CFTs, where the conformal anomaly exists and is parametrized by a single constant c. For higher dimensions, the "twist fields" are nonlocal "line" operators and you need to do some work. Free fields have been studied by Casini and Huerta ( ttps://arxiv.org/abs/0905.2562), and AdS/CFT gives some results from holography. Casini and Huerta have also shown that the entanglement entropy of a circle in a (2+1)-d CFT is equal to the Euclidean free energy on the sphere (see the Hartman lectures you asked about in another recent thread), and there are isolated results for certain CFTs which can be perturbatively accessed using this Calabrese-Cardy replica trick you describe. But general results for general regions in strongly-interacting CFTs are rare.