# Finite element discretization of Reaction-diffusion problem with Dirac source term

I'm writing a code using continuous piecewise linear finite elements on triangular grids to solve the diffusion-reaction problem. the source function f is a Dirac mass at the center. How can i compute the mesh in order to have the dirac in the baricenter of almost one triangle of the mesh? thanks

• What is the domain? You wanna triangulate the domain so that the Dirac impulse function sits on the barycenter of a certain triangle? – Shuhao Cao Nov 27 '13 at 21:08
• the domain is a square of base 2 and I want to triangulate properly becouse if not the solution is 0 – Betelgeuse Nov 27 '13 at 21:54
• What program are you using to do the meshing? It's possible to do this in Triangle. – Daniel Shapero Nov 27 '13 at 22:34
• You can just prescribe the triangle that you want as an interior interface that describe the domain to the mesh generator. Afterwards, you can do 1:4 refine so that the Dirac point is always in the barycenter of a triangular element. FreeFem++, gmsh can do this. – Hui Zhang Nov 28 '13 at 20:45
• At the moment I'm using the function initmesh() of the PDE toolbox – Betelgeuse Nov 29 '13 at 11:17

## 2 Answers

Since you are creating the mesh specifically for this right-hand side in order to have a correct solution, it would make more sense to have the Dirac located on a node instead: On the one hand, this is likely to be easier to realize in the mesh generator, and on the other hand the assembly is easier: Since you are using continuous piecewise linear elements and (presumably) nodal basis functions $\psi_i$, you have for $f=\delta_{x_i}$ (where $x_i$ is one of the nodes of your triangulation) that $$\langle f,\psi_j\rangle_{C^*,C} = \psi_j(x_i) = \begin{cases} 1 &\text{if } i=j,\\ 0 &\text{if } i\neq j, \end{cases}$$ i.e., your load vector is just the vector with $1$ in the $i$th component and $0$ else. (This also works if you have linear combinations of Dirac measures located on the nodes. This also explains what happens when you forget the mass matrix on the right hand side...)

• yes I have the basis functions and I understand your comment, but it's not clear to me that if I locate the Dirac on a node instead of a barycenter, then the solution will be 0 or infinity, or not? – Betelgeuse Nov 29 '13 at 8:29
• For a Dirac on a node, the finite element solution is well-defined (i.e., not "infinity") and not zero (since the linear system arising from your discretization has a unique solution with non-zero right-hand side). – Christian Clason Nov 29 '13 at 8:33

You cannot use the quadrature points to sample the Dirac data. You would either get zero or infinity as the solution. You could choose to use a smooth approximation (e.g. a narrow Gaussian) of the source term if you don't have access to the piece of code performing the integration over elements.

If you can add a "Dirac-loop" to your FEM code, then you can find the entries $f_i$ of the rhs of your linear system $\mathrm{K} \mathbf{u} = \mathbf{f}$, by directly evaluating your shape functions at the dirac points, i.e. $f_i = \int \delta_{\xi_0} N_i = N_i(\xi_0)$. For this, you'd have to determine the local coordinate $\xi_0$ of the dirac load. This can be done by first detecting on which element the dirac load falls (see convex hull algorithms) and then inverting the geometric map of that element either directly (straight elements) or by Newton iteration (curved elements).