Add Dirac to a function

Question

Hi All,

I'm solving a linear evolution equation (ex: heat equation) with a source term that is a dirac in time and space:

du/dt-nabla^2(u)=f(t)dirac(t,x)

In order to do so, i use a finite difference method in time but i need to add some jumps to the solution when it's needed. Let's suppose my time dirac is in t_d.

I use a regular time scheme to solve the homogeneous equation on [0, t_d^-]. To get the solution at time t_d, i need to add a dirac term to my solution (a jump):

u(t_d) = u(t_d^-) + Delta_x f(t)

Is it possible to do so in Fenics ? I'm familiar with the Point-source class, but shall i solve the equation:

y = pointsource(...,magnitude = f(t))

or is there any other solution ?

Thanks by advance,

JayC

RR · Answer 1 · Dec 15, 2016

Hi, I'm using my own Dirac Expressions for Electromagnetic scattering problems. A simple Point source depending on some parameters can be defined as follows:

class X_directed_Source_in_3D(Expression):

    # set your parameters
    # J = Magnitude
    # a = Strech-Factor of Dirac Function
    # origin = [x, y, z] coordiantes of source

    ## Solution 1 (definitely works): give "**df_kwargs" here and "cell=tetrahedron" with calling the Expression instance. Without kwargs, Dolfin prints an "geometric dimension" error during assembly
    ## Solution 2 (referring to question https://fenicsproject.org/qa/3813/expression-subclass-overloading-__init__-requires-element not tested) give "element=None" instead of "**kwargs"

    def __init__(self, J=1.0, a=0.6, origin=[0.0, 0.0, 0.0], **df_kwargs):

        self.J = J
        self.a = a
        self.x = origin[0]
        self.y = origin[1]
        self.z = origin[2]

    def eval(self, values, x):

        r = (x[0] - self.x)**2 + (x[1] - self.y)**2 + (x[2] - self.z)**2
        values[0] = ((J / sqrt(a * pi)) * exp(-r / a**2))
        values[1] = 0.0
        values[2] = 0.0

    def value_shape(self):
        return (3,)

# ... define your J(t) and origin(t)

Solution 1 - call Expression as follows (definitely works):

Dirac_xt = X_directed_Source_in_3D(J=J(t), origin=origin(t), cell=tetrahedron)

Solution 2 - call Expression as follows (not tested, V is your FunctionSpace):

Dirac_xt = X_directed_Source_in_3D(J=J(t), origin=origin(t), element=V.ufl_element())

For other directions, set "values[1]" or "[2]", respectively.

For a 1D FunctionSpace, remove "values[1]" or "[2]" and "def value_shape"
For a 2D FunctionSpace, remove "values[2]" and set "def value_shape" to "return (2,)"

Of course, if you have a 1D or 2D problem, origin has 1, or 2 dim instead of 3 and "r" reduces to

    r = (x[0] - self.x)**2 #or
    r = (x[0] - self.x)**2 + (x[1] - self.y)**2

The parameter surely needs to be adjusted to your problem.
I hope I understood your problem correctly and this solves your issue.

Add Dirac to a function

1 Answer