@tomleslie 

remove(has, EQI__4, 0=0);

@Ramakrishnan Change

v / sqrt(v^+ . v):

to

v / sqrt(v^%T . v):

That should work with older versions of Maple.

@Preben Alsholm Unless I am missing the obvious, there is something wrong with what Maple does with your calculation steps.  To see that, let's give a name to the line that comes after you define L:

K := map(IntegrationTools:-Change,L,cos(2*x)=z);

The entries of L and K should be pairwise equal.  Let's check:

evalf([L[1], K[1]]);  # this gives [0.2898392982,  0.2898392982] which is OK
evalf([L[4], K[4]]);  # this gives [0.2898392982, -0.2898392982] which is WRONG

Vv's assertion that J=0 is still correct, but unless I am mistaken, that does not follow from your calculations.

@brian bovril I spent some time this weekend on developing a procedure which serves to extrapolate a square's boundary data into its interior.  You will find the details in this MaplePrimes post.

The construction there assumes the values of the desired function are given on all four edges of the square.  It does not handle the cases, as in this current thread, where the derivatives are specified on some edges.  I believe that the method I have described in the link noted above can be extended to handle derivative data as well but I haven't attempted that.

@brian bovril I can explain that construction, however it is quite ad hoc and not particularly illuminating.  Here is an alternative solution (of course, the question does not have a unique answer) which is both neater and easier to explain.

Consider the Cartesian coordinate system r−z (r is horizontal, z is vertical) and the associated polar coordinates ρ−θ.  The function f(ρ,θ) = cos θ defined on the square [0,1]×[0,1] is 1 along the bottom edge and 0 along the left edge. That's the same as the function σ₁(r,z) = r/sqrt(r^2+z^2) in the Cartesian coordinates.  Through translation and reflection, we see that the function σ₂(r,z) ≡ σ₁(1-r,1-z), defined on the same square, is 1 along the top edge and 0 along the right edge.

The function σ₃(r,z) ≡ σ₂(r^2,z^2) defined on the larger square [−1,1]×[−1,1] is symmetric with respect to reflection about the coordinate axes, and therefore ∂σ/∂z=0 on the r axis, and  ∂σ/∂r=0 on the z axis.

We observe that the restriction of σ₃(r,z) to the fourth quadrant satisfies the desired boundary conditions.  The only trouble is that the original question places the domain in the first quadrant.  We remedy that through a translation and arrive at the final solution σ(r,z) ≡ σ₃(r,z−1).  Here are the details in Maple:

@ecterrab An excellent job!  As always, I am impressed by your quick response to suggestions made in this forum, and your ever-increasing contributions to Maple.

@rlopez It depends on what one assumes on "r".  Here are two scenarios:

1. We have w = f(b, phi, r), and g(b, phi, r)=0.  The latter may be solved, in principle, for r and thus give r = G(b, phi).  Then w = f(b, phi, G(b,phi) is a function of two independent variables b and phi.
2. We have w = f(b, phi) and g(b,phi)=0.  Then latter may be solved, in principle, for phi and thus give phi = G(b).  Then w = f(b, G(b)) is a function of a single variable b.  Alternatively, we may solve g(b,phi)=0 for b and get b=H(phi).  Then w = f(H(phi), phi) is a function of the single variable phi.

In my calculation I followed the first scenario, while in your calculation you follow the second scenario.  The difference is that in the first scenario r is a variable, while in the second scenario it is a constant. As vv commented, the problem is not clearly formulated.

@ecterrab I am happy to hear that further work is being done on this and I look forward to seeing and using the resulting product.

Thank you everybody for your solutions, comments, and insights.  Now that the problem is solved, I would like to tell you what it is all about.

Consider the standard quadratic equation a*x^2 + b*x + c = 0.  Suppose that the coefficients a, b, c, are selected randomly from a normal (Gaussian) distribution of mean zero and standard deviation sigma. We wish to determine the probability of complex roots.

Any selection a, b, c, of the coefficients corresponds to point (a,b,c) in a three-dimensional Cartesian coordinates system.  The implicitly defined surface b^2 - 4*a*c = 0 divides the space into two subsets, one corresponding to real roots, the other to complex roots.  Here is what the surface looks like:
plots:-implicitplot3d(b^2 - 4*a*c = 0,
    a=-4..4, b=-4..4, c=-4..4,
    style=surface, numpoints=10000,
    orientation=[-65,80,0]);

Aside: Actually that's the surface of a cone, but I am not going to exploit that property here.

It is evident that the region corresponding to complex roots is the union of two subregions
    { (a,b,c) : a > b^2/(4*c), c > 0  }  and  { (a,b,c) : a < b^2/(4*c), c < 0  }.
The two subregions are congruent.  Therefore it suffices to compute the probability of a point lying in one of them and then double the result.

@Mariusz Iwaniuk Thanks for your answer.  Actually, motivated by your success, I played around a bit by changing the orders of the integrations, and thus obtained an exact (non-numeric!) answer which of course agrees with yours:

@DJJerome1976 As to "it never dawned on me that shadebetween( ) accepts variable limits", the documentation may not be without blame.  The help page says:

shadebetween(f1, g1, x=a..b, y=c..d);
where "a, b, c, d are real constants".  (my emphasis)

Furthermore, due to the lack of sufficient documentation, it is not clear why the second of the following two commands fails to "shade in between" in Maple 2018.1:
plots:-shadebetween(0, x^2 + y^2, x=0..1, y=0..x);
plots:-shadebetween(0, x^2 + y^2, y=0..x, x=0..1);

Kitonum's plot3d() variant succeeds in both cases.

@mmcdara You have correctly described the classical Poiseuille flow.  In that context, the radial velocity v is zero and the axial velocity u is parabolic.

The problem originally stated in this thread is not quite the Poiseuille flow—the pipe leaks along its lateral surface as indicated by the boundary condition v=p+a at r=1.  Consequently, the radial velocity is nonzero, and as far as I can see, there is no solution in terms of elementary functions.  A numerical solution should work.

@mmcdara I think the pbar in the original statement is meant to indicate the average pressure.  Thus, the average pressure at z=-1 is zero and at z=+1 is p_a.

I don't know whether these are sufficient boundary conditions to make a well-posed problem.

@Ronan the last command in my worksheet contains the construction
    seq(frame(t), t=0..100, 0.4)

This produces animation frames over the time t=0 to t=100, in steps of 0.4.  You may change those numbers as you wish.

As to the blog that you referred to, I find the interface non-intuitive and therefore am having difficulty navigating through it.  I searched there for your statement regarding the ratio of the angular momentum and kinetic energy but was unable to locate it.  Can you post a direct URL?  The statement should involve more than just the ratio, since that ratio is not a dimensionless quantity.

@vv You are right; it's not difficult to compute moments of inertia of the T-handle exactly through direct integration.  Following your idea, in the attached worksheet I calculate the moments of inertia while accounting for the overlap between the two cylinders.  As expected, the results are quite close to the approximate calculation where the overlap is ignored.

T-handle-flip-extra.mw