# Inactive operator in NDSolve (FEM) produced different solution

GROUPS:
 I was trying to a PDE using FEM. The domain is a cube with Dirichlet boundary conditions imposed( 1 for x=1 face while 0 for the remaining faces). I tried to include the inactive operator to the equation as well, and the results were different. a = {1, 2, 10}; opp = Inactivate[Div[a*Grad[u[x, y, z], {x, y, z}], {x, y, z}], Div | Grad]; uval = NDSolveValue[{opp == 0, u[x, y, -1] == u[x, y, 1] == u[-1, y, z] == u[x, 1, z] == u[x, -1, z] == 0, u[1, y, z] == 1}, u, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, Method -> "FiniteElement"]; opp2 = Div[a*Grad[u[x, y, z], {x, y, z}], {x, y, z}]; uval2 = NDSolveValue[{opp2 == 0, u[x, y, -1] == u[x, y, 1] == u[-1, y, z] == u[x, 1, z] == u[x, -1, z] == 0, u[1, y, z] == 1}, u, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, Method -> "FiniteElement"]; Plot[{uval[x, 0, 0], uval2[x, 0, 0]}, {x, -1, 1}] the first solution actually corresponds to the solution when a={1,1,1}. Is there any way to have the PDE with Inactive operator solved correctly? Since I am going to replace the variable a by a discontinuous function, evaluating its derivative before passing to NDSolveValue is not desirable (see section for inactive operator in http://reference.wolfram.com/language/tutorial/NDSolvePDE.html ). Attachments: