Position of discontinuous coefficient influences the solution of PDE

Here is an explanation of what happens. Let's setup the problem once more.

Lmid = 1; L = 2; tend = 1;
m[x_] = If[x < Lmid, 1, 2];
(*m[x_]=2;*)
eq1 = m[x] D[u[x, t], t] == D[u[x, t], x, x];
eq2 = D[u[x, t], t] == D[u[x, t], x, x]/m[x];
opts = Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}};
ndsolve[eq_] := 
  NDSolveValue[{eq, u[x, 0] == Exp[x]}, u, {x, 0, L}, {t, 0, tend}, 
   opts];

Equation 1 and 2 are mathematically the same, however, when we evaluate them we get different results as shown here:

sol1 = ndsolve[eq1];
Plot[sol1[x, tend], {x, 0, L}]

sol2 = ndsolve[eq2];
Plot[sol2[x, tend], {x, 0, L}]

What happens? Let's look at how the PDE gets parsed.

ClearAll[getEquations]
getEquations[eq_] := Block[{temp},
  temp = NDSolve`ProcessEquations[{eq, u[x, 0] == Exp[x]}, 
     u, {x, 0, L}, {t, 0, tend}, opts][[1]];
  temp = temp["FiniteElementData"];
  temp = temp["PDECoefficientData"];
  (# -> temp[#]) & /@ {"DampingCoefficients", "DiffusionCoefficients",
     "ConvectionCoefficients"}
  ]

getEquations[eq1]
{"DampingCoefficients" -> {{If[x < 1, 1, 2]}}, 
 "DiffusionCoefficients" -> {{{{-1}}}}, 
 "ConvectionCoefficients" -> {{{{0}}}}}

This looks good.

getEquations[eq2]
{"DampingCoefficients" -> {{1}}, 
 "DiffusionCoefficients" -> {{{{-(1/If[x < 1, 1, 2])}}}}, 
 "ConvectionCoefficients" -> {{{{-(If[x < 1, 0, 0]/
       If[x < 1, 1, 2]^2)}}}}}

For the second eqn. we get a convection coefficient term. Why is that? The key is to understand that the FEM can only solve this type equation:

$d\frac{\partial }{\partial t}u+\nabla \cdot (-c \nabla u-\alpha u+\gamma ) +\beta \cdot \nabla u+ a u -f=0$

Note, that there is no coefficient in front of the $\nabla \cdot (-c \nabla u-\alpha u+\gamma)$ term. To get things like $h(x) \nabla \cdot (-c \nabla u-\alpha u+\gamma)$ to work, $c$ is set to $h$ and $\beta$ is adjusted to get rid of the derivative caused by $\nabla \cdot (-c \nabla u)$

Here is an example:

c = h[x];
β = -Div[{{h[x]}}, {x}];
Div[{{c}}.Grad[u[x], {x}], {x}] + β.Grad[u[x], {x}]
(* h[x]*Derivative[2][u][x] *)

In the case at hand that leads to:

Div[{{1/m[x]}}.Grad[u[x], {x}], {x}] - 
  Div[{{1/m[x]}}, {x}] // Simplify

(* {Piecewise[{{Derivative[2][u][x]/2, x >= 1}}, Derivative[2][u][x]]} *)

But that is the same as specifying:

 eq3 = D[u[x, t], t] == 
   Inactive[
     Div][{{1/If[x < 1, 1, 2]}}.Inactive[Grad][u[x, t], {x}], {x}];

sol3 = ndsolve[eq3];
(* Plot[sol2[x, tend] - sol3[x, tend], {x, 0, L}] *)

I have checked that flexPDE (another FEM tool) gives exactly the same solutions in all three cases. So this issue is not uncommon. In principal a message could be generated but how would one detect when to trigger that message? If you have suggestions about this, let me know in the comments. I think it were also good to add this example to the documentation - if there are no objections. I hope this clarifies the unexpected behavior a bit.