Implement point source in the variational formulation

variational formulation
1

Dear FEniCS community,

I’m having a hard time trying to implement a point load on a 3D problem in 1D mesh. I have read old fenics post, i.e., https://fenicsproject.org/qa/14258/how-to-apply-pointsource-to-linearvariationalproblem/ or https://fenicsproject.org/qa/7941/applying-dirac-function-using-fenics-pointsource-function/ , with no luck.

The sample code is as simple as the following:

DOF = 50

# Create mesh and define function space
L=4
mesh = IntervalMesh(DOF, 0, L)

degree = 1
U = VectorFunctionSpace(mesh, 'CG', degree, dim=3)

u = TrialFunction(U)

v = TestFunction(U)

# Dirichlet on left boundary
def l_boundary(x, on_boundary):
    return on_boundary and near(x[0], 0.0)

zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(U, zero, l_boundary)

# Variational formulation
A = inner(grad(u), grad(v))*dx

# I want the following load f to be a point load on the endpoint L=4
L = inner(f, v)*dx 

# Solution to the problem
u_h = Function(U)
solve(A == L,u_h,bc)

From the post related above, a suggested solution was that an approximation of the Dirac delta function can do the work, but I have no idea on how to do this on FEniCS. I tried the solution from Mirok but get an Maximum recursion depth exceeded error.

Thanks in advance for any help!.

2

The Expression approach should be modified as follows

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

    def value_shape(self): return (3, )
3 Likes
3

Dear Mirok,

Thank your for your solution. The code runs succesfully if I set delta=Delta(1E-4) and then

L=inner(delta,v)*dx

But I don’t know how this works, i.e, ¿how to set the magnitude and location of the source given that location is x=4 and source is f=(0,0,4)?

4

MiroK has outlined the code for a mollified Dirac delta function. See here.

For the precise formulation of a point source in your variational formulation, and FEniCS code, use the PointSource class documented here.

5

Dear nate,

I missed to point out in my question that I’m trying NOT to use the PointSource class in order to avoid the assemble() and apply() part. That will mess up my code…

However, thanks for the documentations.

6

Changing location mounts to

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

delta = Delta(eps=1E-4, x0=np.array([0.5, 0, 0]), degree=5)

You can control magnitude by scaling, i.e. L = inner(Constant(100)*delta, v)*dx

4 Likes
7

Dear MiroK,

Thank you for your valuable help. The code works successfully. For a 1D mesh, we can call Delta as:

delta=Delta(eps=1E-4, x0=0.5, degree=5)
9

Hello @MiroK Mirok, what is the magnitude of the force and location in this code?
I want to add a point load at the free end of a cantilever, is this code useful for that?
Thanks.

10

Hi Mirok, I have a question to ask you. When I run this code, it returns this error. Is this a problem with my “fenics_adjoint”, or is it something else?

fig1
fig11271×193 82.3 KB

11

Using user-expressions in dolfin_adjoint requires Extra attention, see

You also need to actually post your code, encapsulated with 3x`, as

```python
from fenics import *
from fenics_adjoint import *
# add more code here
```

and make sure that the code is a minimal reproducible example.

1 Like
12

Thank you very much for your help, I will read these and find my mistakes.

13

Hi, this seems to be an issue more related to dolfin-adjoint. See here https://bitbucket.org/dolfin-adjoint/pyadjoint/issues/122/update-expression-values-with-solution

14

Dear MiroK, I think I have a new question. The example is applying a load to a single node, how can I load multiple nodes at different locations at the same time? Thank you very much for your answer.

15

Hi, one way building on the previous code is

deltas = [Delta(eps=1E-4, x0=np.array([0.5, 0, 0]), degree=5),
          Delta(eps=1E-4, x0=np.array([0.6, 0, 0]), degree=5),
          Delta(eps=1E-4, x0=np.array([0.7, 0, 0]), degree=5)]

L = sum(inner(Constant(100)*delta, v)*dx for delta in deltas)
2 Likes