How to solve wave equation with finite elements when material properties vary continuously over a region

Question

I am solving the one-dimensional wave equation over regions where the bulk modulus (and thus the wave speed) vary continuously over a region. The current version seems to assume that the material properties are constant over an FEM element. Some literature suggests that isoparametric elements can help model regions where the material properties vary continuously. Can you suggest how to model regions where the material properties vary in a continuous fashion?

In the wave equation code I am using, kappa and rho can vary across an element. The code is as follows:

eqn = 1/κ[x] D[u[t, x], {t, 2}] + 
Div[-ρ[x]*Grad[u[t, x], {x}], {x}] == 
1/κ[x]*10*Exp[-50 (x^2)]*Sin[2 π f t] + 
NeumannValue[0, x == 0] + NeumannValue[-Derivative[1, 0][u][t, x], x == 
xMax]; 
ic = {u[0,x] == 0, Derivative[1, 0][u][0, x] == 0};

Answer 1

Here is a wave equation with a $x+1$ as a factor in the wave equation.

\[CapitalOmega] = 
  RegionDifference[
   RegionDifference[Rectangle[{0, 0}, {2, 1}], 
    Rectangle[{9/10, 0}, {11/10, 4/10}]], 
   Rectangle[{9/10, 6/10}, {11/10, 1}]];
sol = NDSolveValue[{D[u[t, x, y], {t, 2}] - 
      Inactive[Div][(x + 1) Inactive[Grad][u[t, x, y], {x, y}], {x, 
        y}] == 0, DirichletCondition[u[t, x, y] == 0, True], 
    u[0, x, y] == 2*Exp[-125 ((x - 0.25)^2 + (y - 0.5)^2)], 
    Derivative[1, 0, 0][u][0, x, y] == 0}, 
   u, {t, 0, 2}, {x, y} \[Element] \[CapitalOmega]];

ListAnimate[
 Table[Rasterize[
   Plot3D[sol[t, x, y], {x, y} \[Element] \[CapitalOmega], 
    PlotRange -> {-0.75, 2}, AspectRatio -> Automatic]], {t, 0, 2, 
   1/25}], SaveDefinitions -> True]

I would be interested to know what lead you to the conclusion that "...current version seems to assume that the material properties are constant over an FEM element." This is not correct and if this appears or is suggested in the documentation it needs to be fixed. But for that I'd need to know what lead you to that conclusion.