Context
-------

I would like to compute the torque that a  (thin) disc applies onto a ring. I.e. I try to understand what is the impact of this outer ring on
the inner disc in the simulation below.

![Mathematica graphics](http://i.stack.imgur.com/DZtTP.png)

For this I would like to *compute the gravitational potentia*l generated by a (razor) thin disc and a ring.
So the abstraction is the following (seen from 9 different angles)

![Mathematica graphics](http://i.stack.imgur.com/jWVoS.png)

Once I know how to compute the potential it becomes straightforward to compute the torque one feature applies on the other.
I want to use FEM here for flexibility when I will account for a more realistic abstraction of the problem.


Attempt
-------

I have defined a domain

    dom = ImplicitRegion[0 <= x <= 1 && -1 <= y <= 1, {x, y}];
 
and the Laplace operator
   
    op = -Laplacian[u[x, y], {x, theta, y}, "Cylindrical"];

I impose the edge condition that the potential should be 1 on the disc
   

    edge = {DirichletCondition[u[x, y] == 1, 0 <= x <= 1/2 && y == 0};

When I solve 

    uD = NDSolveValue[{op == 0, edge},  u, {x, y} \[Element] dom]

I get  

    StreamPlot[-{D[#, x], D[#, y]} &@uD[x, y] // Evaluate, {x, y} \[Element] dom]
![Mathematica graphics](http://i.stack.imgur.com/Qw5o3.png)


> **Problems**: 

> (i) The outer box imposes a (box like) symmetry which is not in the sought solution 

> (ii) Strangely enough The code fails if I use the ring-like boundary condition instead:

    {DirichletCondition[u[x, y] == 1, 1/2 <= x <= 3/4 && y == 0]};


Question
--------

> How to compute the gravitational potential created by a disc (and a ring) using FEM in NDSolve?

In a broader sense I think I am asking how can FEM methods deal with PDEs with boundaries at infinity?
*I am guessing that one strategy might be to move the boundary sufficiently far away and increasing sampling within the inner region?*

Note that my attempt above is imposing fixed potential on the disk not fixed density. I am not sure this is important or not, but ideally (to compare to the analytical solution below) fixing density would be better.


**PostScriptum**

I have found this (nice!) [blog][1] which provides me with an analytic solution 
as follows

    PhiDiskData = 
      WolframAlpha[
        "electric potential of a charged disk", {{"Result", 1}, 
         "Input"}] // ReleaseHold;
    PhiDisk = PhiDiskData /.  QuantityVariable[a_, _] -> a /. { 
        Q -> Pi R^2, "ElectricConstant" -> 1};
    Phi = PhiDisk /. { x -> r Cos[Theta], y -> r Sin[Theta]} //
       Simplify[#, Assumptions -> {r > 0}] &

![Mathematica graphics](http://i.stack.imgur.com/N0l1d.png)

    Clear[fD]; fD = 
     FullSimplify[-D[PhiDisk, {{x, y, z}}]  /. x^2 + y^2 -> r^2, 
       Element[z, Reals] && r > 0]  /. {x -> r Cos[Theta], y -> r Sin[Theta]};
       fD = -{Sqrt[fD[[1]]^2 + fD[[2]]^2] // FullSimplify, fD[[3]]};

![Mathematica graphics](http://i.stack.imgur.com/26Cjq.png)

So that 

    phiN = (Phi /. { Theta -> 0, R -> 1/2}); pl1 = 
     ContourPlot[Evaluate[phiN], {r, 0, 2}, {z, -2, 2}, Exclusions -> {}, 
      Contours -> 15,ColorFunction -> (ColorData["RedBlueTones"][1 - #] &),             
      Epilog -> {Thickness[0.02], Line[{{0, 0}, { 1/2, 0}}]},  
      FrameLabel -> {r, z}, PlotRange -> All, AspectRatio -> 2];
    pl2 = StreamPlot[(fD /. R -> 1/2), {r, 0, 2}, {z, -2, 2},
       AspectRatio -> 2, StreamStyle -> White];
    pl3 = Show[pl1, pl2, PlotRange -> {{0, 1.5}, {-0.5, 1}}, 
      AspectRatio -> 1]

yields

![Mathematica graphics](http://i.stack.imgur.com/bUX03.png)

> So my question amounts to finding this solution numerically.


Note that the analytic solution works nicely for rings as well (if defined as the difference between two discs.)

    phiN = (Phi /. { Theta -> 0, 
         R -> 1}) - (Phi /. { Theta -> 0, R -> 1/2}); pl1 = 
     ContourPlot[Evaluate[phiN], {r, 0, 2}, {z, -2, 2}, Exclusions -> {}, 
      Contours -> 15,ColorFunction -> (ColorData["RedBlueTones"][1 - #] &),
    Epilog -> {Thickness[0.02], Line[{{1/2, 0}, { 1, 0}}]},  
      FrameLabel -> {r, z}, PlotRange -> All, AspectRatio -> 2];
    pl2 = StreamPlot[(fD /. R -> 1) - (fD /. R -> 1/2), {r, 0, 2}, {z, -2,
         2},AspectRatio -> 2, StreamStyle -> White];
    pl4 = Show[pl1, pl2, PlotRange -> {{0, 1.5}, {-0.5, 1}}, 
      AspectRatio -> 1]

![Mathematica graphics](http://i.stack.imgur.com/RCWcD.png)


  [1]: http://blog.wolfram.com/2012/09/27/3d-charges-and-configurations-with-sharp-edges/