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

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}] +
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};

• I don't understand your question. Have you tried to define example function kappa and rho and to submit the system to NDSolve? What is the issue? – Henrik Schumacher Aug 5 at 1:45
• "The current version seems to assume that the material properties are constant over an FEM element." Might be true for element order 1. For order 2, probably a nonexact quadrature rule is applied, but I doubt that the coefficients are assumed to be constant. Anyways, with sufficiently small elements and suffciently smooth material coefficients, that should not matter at all. – Henrik Schumacher Aug 5 at 1:49
• Thanks, I use order 2 elements. Some literature suggests that accuracy is improved by using isoparametric elements.This shows up in functionally graded materials (FGM). My current results seem OK, using a fine grid in the region where the coefficients are rapidly varying. – Klandgren Aug 6 at 13:34

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.