Dirac point source class implementation for Poisson equation

Legacy DOLFIN
1

Hello,

I have been trying to adapt the proposed Delta (source point) class by “MiroK” in this post Delta variational formulation, but I keep getting a shape mismatch error. My code is :

class Delta(UserExpression):
    def __init__(self, eps, x0, **kwargs):
        self.eps = eps
        self.x0 = x0
        UserExpression.__init__(self, **kwargs)
    def eval(self, values, x):
        eps = self.eps
        values[0] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
        values[1] = 0
        values[2] = 0

    def value_shape(self): return (3, )

from mshr import *
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

resolution = 15
domain_geometry = Circle(dolfin.Point(0.5, 0.5), 0.5+0.01)
mesh = generate_mesh(domain_geometry, resolution)
V = FunctionSpace(mesh, "Lagrange", 1)

def boundary(x, on_boundary):
    return on_boundary

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
delta = Delta(eps=1E-4, x0=np.array([0.5, 0, 0]), degree=5)
f = delta
magnitude = Constant(100)
g = Constant(0.0)
a = inner(grad(u), grad(v))*dx
L = inner(magnitude*f, v)*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

plot(u)
plt.colorbar(plot(u))

The class was initialy done for a 3 dimensions problem. Since i’m working with a 2D problem, how can I modify the class, or maybe adapt my code to it.

The error I keep getting is on line : L = inner(magnitude*f, v)*dx
Saying : shapes do not match.

Thank you.

2

Hi

The error comes from the fact that your are solving a 2D problem and using a scalar source. You have to redefine the Class to be scalar-type. The following works and impose a point source at (0.5,0.5)

from mshr import *
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np


class Delta(UserExpression):
    def __init__(self, eps, x0, **kwargs):
        self.eps = eps
        self.x0 = x0
        UserExpression.__init__(self, **kwargs)

    def eval(self, values, x):
        eps = self.eps
        values[0] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)

    def value_shape(self): return ()


resolution = 15
domain_geometry = Circle(dolfin.Point(0.5, 0.5), 0.5+0.01)
mesh = generate_mesh(domain_geometry, resolution)
V = FunctionSpace(mesh, "Lagrange", 1)


def boundary(x, on_boundary):
    return on_boundary


# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
delta = Delta(eps=1E-4, x0=np.array([0.5, 0.5]), degree=5)
f = delta
magnitude = Constant(100)
g = Constant(0.0)
a = inner(grad(u), grad(v))*dx
L = inner(magnitude*f, v)*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

Cheers

2 Likes
3

Thank you ! It works :slight_smile: !