# 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

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...)
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).