Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It only takes a minute to sign up.
Sign up to join this community
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.
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.
