# Computing element stiffness matrices with variable coefficients

I am trying to implement a simple FEM approach, using p1 triangular elements, for solving the diffusion equation with variable nodal diffusivities and I was wondering how to incorporate the variable nodal diffusivities when computing the element stiffness matrices.

I would really appreciate if someone could point me in the right direction. Many thanks in advance.

I assume you're trying to solve an equation that looks like:

\begin{align} -\nabla \cdot (a(x)\nabla{u}) = f, \end{align}

for $x$ in some domain $\Omega$, although the same approach would be fine (for a residual evaluation, anyway) if $a$ were also a function of $u$.

The stiffness matrix will take the form

\begin{align} A_{ij} = \int_{\Omega}a(x)\nabla\varphi_{j}(x)\cdot\nabla\varphi_{i}(x)\,\mathrm{d}x. \end{align}

For $P_{1}$ triangular elements, $\nabla\varphi_{j}$ will be constant over each element $K$, for each $j$. Letting $\mathcal{K}$ be the mesh, the elements of the stiffness matrix become

\begin{align} A_{ij} = \sum_{K \in \mathcal{K}}\nabla\varphi_{j}^{K}\cdot\nabla\varphi_{i}^{K}\int_{K}a(x)\,\mathrm{d}x, \end{align}

where all I've done is factor out the constant gradients over each element (these gradients will differ from element to element, which is why they are also indexed over $K$). If $K$ has vertices $N_{1}$, $N_{2}$, and $N_{3}$, then

\begin{align} \int_{K}a(x)\,\mathrm{d}x = a\left(\frac{N_{1} + N_{2} + N_{3}}{3}\right)|K| + O(h^{3}), \end{align}

where $h$ is the mesh "size", and $|K|$ is the area of element $K$. This approximation is also called the center of gravity rule.

In practice, you'd take the local stiffness matrix for each element (assuming a diffusivity of one), and multiply it by the variable diffusivity $a$ evaluated at the center of gravity of the element.

• Hi Geoff, many thanks for taking the for an excellent reply. If I may pose a following question: is there a case when it is necessary to expand these variable diffusivities in terms of the nodal basis functions? Many thanks again. – semper May 21 '14 at 11:28

If you are using numerical integration you can simply replace the integral $$\int\limits_{\Omega_{el}} a(x) \nabla \varphi_i(x) \nabla \varphi_j(x) dx = \sum_{x_k \in \Omega_{el}} a(x_k) \nabla \varphi_i(x_k) \nabla \varphi_j(x_k) \enspace .$$

So, you have to know your diffusivity at each Gauss point. This is easily know if you have an analytic function. If you know the data just in your nodes you can use the interpolators yo compute these values. In the case of linear elements is simpler to do what Geoff suggested.