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: