# Problem with Neumann condition in quarter disc

Bug introduced in 10.0 and fixed in 10.2

So I'm following the available examples in version 10 for FEM, The plane stress operator is shown as this

op0 = {Inactive[
Div][({{0, -((Y ν)/(1 - ν^2))}, {-((Y (1 - ν))/(
v[x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-(Y/(1 - ν^2)),
0}, {0, -((Y (1 - ν))/(2 (1 - ν^2)))}}.Inactive[
Grad][u[x, y], {x, y}]), {x, y}],
Inactive[
Div][({{0, -((Y (1 - ν))/(2 (1 - ν^2)))}, {-((
u[x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-((Y (1 - ν))/(2 (1 - ν^2))),
v[x, y], {x, y}]), {x, y}]}


 op = op0 /. {Y -> 10^3, ν -> 33/100};


My goal is to model a disc with pressure in the outer surface. I placed a boundary condition that reduced the computational domain to the upper quarter of the plane:

Subscript[Γ, D] = {
DirichletCondition[{u[x, y] == 0.}, x == 0],
DirichletCondition[{ v[x, y] == 0.0}, y == 0]}

ℛ =
ParametricRegion[{{s, t},
Sqrt[s^2 + t^2] <= 10.0 && 0 <= t <= 10.0 && 0 <= s <= 10.0}, {s,
t}];

dd = DiscretizeRegion[ℛ]


At this point I get as expected, a mesh in the upper quarter plane of the disc of radius 10.0

Now I try to solve the system applying some Neumann condition in the outer radius, but apparently something about it is failing it:

{uif, vif} =
NDSolveValue[{op == {
NeumannValue[-1.3*(x/10.0), Sqrt[x^2 + y^2] >= 9.9],
NeumannValue[-1.3*(y/10.0), Sqrt[x^2 + y^2] >= 9.9]},
Subscript[Γ, D]}, {u, v}, Element[{x, y}, dd]];


I get

Power::infy: Infinite expression 1/Sqrt[0.] encountered. >>
Infinity::indet:Indeterminate expression 0. ComplexInfinity encountered. >>
CompiledFunction::cfta: Argument {Boole[-9.9>=Indeterminate],0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,<<14>>} at position 1 should be a rank 1 tensor of machine-size integers. >>
CompiledFunction::cfta: Argument {Boole[-9.9>=Indeterminate],0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,<<14>>} at position 1 should be a rank 1 tensor of machine-size integers. >>

• Yes, this is a bug and it's filed. Luckily you found a good workaround. One small other thing: If you want to use a mesh as input to NDSolve I'd use ToElementMesh[R] to create the mesh. That mesh will be second order (more accurate) while DiscretizeRegion produces a first order mesh. Nov 5, 2014 at 9:49

Apparently NeumannValue doesn't like the square root, replacing it with the squared expression works correctly:

NeumannValue[1000.3* (x/10.0), x^2 + y^2 >= 99.9]
`