Gravitational potential created by a thin disc using FEM and NDSolve

The method outlined by Henrik Schumacher can also be using while explicitly relying on cylindrical symmetry (when applicable), so as to reach better accuracy.

Laplace equation with a given boundary

    Needs["NDSolve`FEM`"];
dz = 1/32; reg0 = Rectangle[{{1/2, -dz}, {1, dz}}];
reg1 = RegionDifference[Disk[{0, 0}, 4], 
   Rectangle[{1/2, -dz}, {1, dz}]];
reg = RegionIntersection[reg1, ImplicitRegion[x > 0, {x, y}]];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.01, 
   "MaxBoundaryCellMeasure" -> 0.025];
Show[mesh["Wireframe"], Frame -> True, PlotRange -> All]

Then the solution reads

sol = NDSolveValue[{-Laplacian[u[x, z], {x, y, z}, "Cylindrical"] == 
      NeumannValue[1, 1/2^2 <= x^2 <= 1 && -dz <= z <= dz], 
     DirichletCondition[u[x, z] == 0, x^2 + z^2 >= 2^2]} // Evaluate, 
   u, {x, z} \[Element] mesh];
Off[InterpolatingFunction::dmval]

We can plot it as

plr = Region@ImplicitRegion[1/2^2 <= x^2 <= 1 && -dz < z < dz, {x, z}];
Clear[sol1]; sol1[x_, z_] = sol[Abs[x], z];
plc = ContourPlot[sol1[x, z], {x, -2, 2}, {z, -2, 2}, Contours -> 20, 
   AspectRatio -> 1, ColorFunction -> "Heat"];
pls = Show[
   StreamPlot[{D[sol[Abs[x], z], x], D[sol[Abs[x], z], z]} /. 
      Abs'[x] -> If[x >= 0, 1, -1] // Evaluate, {x, -2, 2}, {z, -2, 
     2}, StreamStyle -> Gray], plr];
Show[plc, pls]    

Poisson equation with a given surface density

Similarly, adapting the 3D case to the cylindrical symmetry

reg = Disk[{0, 0}, 2];
reg = RegionIntersection[reg, ImplicitRegion[x > 0, {x, y}]];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.001];
\[Rho] = Function[{x, z}, 
   Boole[1/2^2 <= x^2 <= 9/4 && -dz <= z <= dz]];
sol = NDSolveValue[{-Laplacian[u[x, z], {x, \[Theta], z}, 
        "Cylindrical"] == \[Rho][x, z] // Evaluate, 
    DirichletCondition[u[x, z] == 0, x^2 + z^2 >= 2^2]}, 
   u, {x, z} \[Element] mesh];
Off[InterpolatingFunction::dmval]
plc = ContourPlot[sol[x, z], {x, -2, 2}, {z, -2, 2}, Contours -> 20, 
   AspectRatio -> 1, ColorFunction -> "Heat"];
pls = Show[
   StreamPlot[{D[sol[Abs[x], z], x], D[sol[Abs[x], z], z]} /. 
      Abs'[x] -> If[x >= 0, 1, -1] // Evaluate, {x, -2, 2}, {z, -2, 
     2}, StreamStyle -> Gray], plr];
Show[plc, pls]

We can check the accuracy of the solution as follows First compute the Laplacian of the solution

f[x_, z_] = -Laplacian[sol[x, z], {x, \[Theta], z},];

then

 {Plot[{f[x, 0.], \[Rho][x, 0.]}, {x, 0, 2}],
  Plot[{f[3/4, x], \[Rho][3/4, x]}, {x, -1, 1}, PlotRange -> All]}

and

{Plot[{f[x, 0.] - \[Rho][x, 0.]}, {x, 0, 2}],
 Plot[{f[3/4, x] - \[Rho][3/4, x]}, {x, -1, 1}, PlotRange -> All]}

I am not too familiar with gravitational potentials but here are some suggestions:

1. I would gauge them to be 0 at infinity and to satisfy Neumann conditions on solid matter.

2. Moreover when specifying boundary conditions with inequalities, you should avoid < and > because the boundary conditions will quite likely not be set correctly, if the boundary is given by the these inequalties with < and > replaced by =.

Afterwards, I gain this, which seems quite plausible to me:

Needs["NDSolve`FEM`"];
dz = 1/16; reg0 = 
 RegionDifference[Cylinder[{{0, 0, -dz}, {0, 0, dz}}, 1], 
  Cylinder[{{0, 0, -dz}, {0, 0, dz}}, 1/2]];
reg = RegionDifference[Ball[{0, 0, 0}, 2], reg0];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.01, 
   "MaxBoundaryCellMeasure" -> 0.025];
sol = NDSolveValue[{
    Laplacian[u[x, y, z], {x, y, z}] == 
     NeumannValue[1, 1/2^2 <= x^2 + y^2 <= 1 && -dz <= z <= dz],
    DirichletCondition[u[x, y, z] == 0, x^2 + y^2 + z^2 >= 2^2]
    }, u, {x, y, z} \[Element] mesh
   ];
Off[InterpolatingFunction::dmval]
GraphicsRow[
 {
  ContourPlot[sol[x, 0, z], {x, -2, 2}, {z, -2, 2}, AspectRatio -> 1, 
   Contours -> 20], 
  ContourPlot[sol[x, y, 0], {x, -2, 2}, {y, -2, 2}, AspectRatio -> 1, 
   Contours -> 25]
  }
 ]

The stream lines of the gravitational force:

plr = DiscretizeRegion[
   ImplicitRegion[1/2^2 <= x^2 <= 1 && -dz < z < dz, {x, z}], 
   PlotTheme -> "Minimal"];
g = Show[
  StreamPlot[-{D[sol[x, y, z], x], D[sol[x, y, z], z]} /. y -> 0 // 
    Evaluate, {x, -1, 1}, {z, -1, 1}, StreamPoints -> Fine
   ],
  plr]

Finally, you could try to apply transformations of following type to your PDE (and the thin disc region) in order to map the sphere at infinity to a sphere of finite radius R:

Φ = X \[Function] X /(R - Sqrt[X.X]);
Ψ = Y \[Function] R Y/(1 + Sqrt[Y.Y]);

Note that you have to transform the Laplacian and also the Neumann conditions; that might be not too convenient...

It might be even more realistic to model the ring as a mass density ρ (for example constant on the ring) instead of cutting the rings out of the computational domain. So you would have to put it on the right hand side of the Laplacian instead of the NeumannValue. This could look like this:

reg = Ball[{0, 0, 0}, 2];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.001];
ρ = Function[{x, y, z}, Boole[1/2^2 <= x^2 + y^2 <= 1 && -dz <= z <= dz]];
sol = NDSolveValue[{
    Laplacian[u[x, y, z], {x, y, z}] == ρ[x, y, z],
    DirichletCondition[u[x, y, z] == 0, x^2 + y^2 + z^2 >= 2^2]
    }, u, {x, y, z} ∈ mesh
   ];
Off[InterpolatingFunction::dmval]
g = GraphicsRow[
  {
   ContourPlot[sol[x, 0, z], {x, -2, 2}, {z, -2, 2}, AspectRatio -> 1,
     Contours -> 20], 
   ContourPlot[sol[x, y, 0], {x, -2, 2}, {y, -2, 2}, AspectRatio -> 1,
     Contours -> 25]
   },
  ImageSize -> Large
  ]

The stream lines of the gravitational force are not perpendicular to the boundary of the disk anymore, but that cannot be expected with a mass density:

plr = DiscretizeRegion[
   ImplicitRegion[1/2^2 <= x^2 <= 1 && -dz < z < dz, {x, z}], 
   PlotTheme -> "Minimal"];
g = Show[
  StreamPlot[-{D[sol[x, y, z], x], D[sol[x, y, z], z]} /. y -> 0 // 
    Evaluate, {x, -2, 2}, {z, -2, 2}, StreamPoints -> Fine
   ],
  plr]

This is in quite a good accordance with the previous result. Note that the used mesh is coarser than the one before, in particular in the region around the ring. One has to fine tune it, e.g. with a suitable choice of a MeshRefinementFunction.

I'm very sorry to spoil this party but solving PDE in the whole volume is a bad idea if all that you want is a force on one body. Use green function $V(\vec{x})=-\underset{\mathbb{R}^{3}}{\int}\frac{G}{|\vec{x}-\vec{r}|}\rho(\vec{r})dr$ to get potential and force exactly in points of the second body and calculate all you need.