Solving Poisson equation with variable coefficient with Finite element method

Question

I need to solve the poisson equation with varibale coefficient $\nabla\cdot(k\nabla u)=f$ with the finite element method (FEM). I really don't know how to obtain the weak formulation of the problem.

Answer 1

The weak form of the problem is found by taking the $L^2$ Inner product with a set of test functions we usually call $V = \{v\in L^2(\Omega): \nabla v\in L^2(\Omega)\}$. This notation defines a Sobolev Space $V$ which is useful for proofs regarding convergence. The important thing to note is that it requires that over the domain $\Omega$, all $v$ are $L^2$ integrable and all elements of the gradient of $v$ are $L^2$ integrable.

The $L^2$ inner product of two 1D functions $u$ and $v$ is defined as: $$ \langle u, v\rangle = \int_a^bu(x)v(x)dx $$ A function $v$ is $L^2$ integrable if over the entire domain $\langle v, v\rangle$ is bounded.

For multidimensional functions $u, v : \mathbb{R}^3\mapsto\mathbb{R}$ the $L^2$ inner product is the volume integral over the domain which I'll notate as: $$ \langle u, v\rangle = \int_\Omega u v \,d\Omega $$

So now to construct the weak form our first step is to take the $L^2$ inner product of the original problem with the test function $v$ and say we have to: $$ \text{find $u$ such that:}\\ \int_\Omega v \nabla\cdot(k\nabla u)d\Omega = \int_\Omega vfd\Omega\quad \forall v\in V $$ Typically we want $u$ and $v$ to both be in the same Sobolev space definition, and right now this form has 2nd derivatives of u in the $L^2$ inner product which isn't garuanteed to be bounded by our definition. So we transform the left hand side using integration by parts. Depending on how you want to view it you can also say that we are using Green's First Identity or the Divergence Theorem with the Product Pule. $$ \int_\Omega (k\nabla u)\cdot\nabla v \,d\Omega = \int_\Gamma v(k\nabla u)\cdot \mathbf{n} \,d\Gamma - \int_\Omega v\nabla\cdot(k\nabla u)\,d\Omega $$

Where $\Gamma$ is the boundary of the domain $\Omega$ and $\mathbf{n}$ is the normal of the boundary. I'll let you do the substitution of the left hand side of the weak form to get the final weak form.

One Note: Depending on the type of Finite Element Method you are using the boundary integral can take on slightly different meanings. If you are using continuous finite elements (look up Continuous Galerkin Methods) then the boundary integral will set by your boundary conditions. This is by no means the only resource but one reference for the treatment of the boundary conditions for continuous finite elements in chapter 9 of Numerical Solution of Differential Equations, Introduction to Finite Difference and Finite Element Methods by Zhilin Li, Zhonghua Qiao and Tao Tang. For discontinuous finite element methods you have to consider the boundaries between elements and usually involves a numerical flux approximation.