# Implement indicator function in the weak formulation

340

views

0

Hello everyone,

I have a term in my diffusion pde which represents the area where my function h is non-zero. I.E. ind(u) = 1 if u>=0 and ind(u) = 0 otherwise.

I have tried a couple of things and cannot figure out how to fit this into FEniCS. See the following:

def ind(u_n):

if u_n >=0:

return 1

else:

return 0

F = u*v*dx + dt*q(u)*dot(grad(u), grad(v))*dx - u_n*v*dx + ind(u_n)*(u_n)*v*dx

doesn't work, neither does

F = u*v*dx + dt*q(u)*dot(grad(u), grad(v))*dx - u_n*v*dx + E*u_n[u_n>0]*v*dx

Can anyone help me? Full code attached at the end.

Thanks, Zhen.

from fenics import *

import time

import scipy

def q(u):

"return nonlinear coefficient"

return u**3

T = 2.0 # final time

num_steps = 20 # number of time steps

dt = T / num_steps # time step size

E = 0.1 # evaporation constant

# Create mesh and define function space

nx = ny = 30

mesh = RectangleMesh(Point(-3, -3), Point(3, 3), nx, ny)

V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition

def boundary(x, on_boundary):

return on_boundary

bc = DirichletBC(V, Constant(0), boundary)

# Define initial value

u_0 = Expression('x[0]>=-1 && x[0]<=1 && x[1] >=-1 && x[1] <= 1 ? exp(-3*a*pow(x[0], 2) - a*pow(x[1], 2)):0',

degree=2, a=5)

#u_0 = Constant(0)

u_n = interpolate(u_0, V)

# Define variational problem

u = Function(V)

v = TestFunction(V)

f = Constant(0)

#note whether there is evaporation or not

def ind(u_n):

if u_n >=0:

return 1

else:

return 0

F = u*v*dx + dt*q(u)*dot(grad(u), grad(v))*dx - u_n*v*dx + E*u_n[u_n>0]*v*dx

# Create VTK file for saving solution

vtkfile = File('two_d_spreading/solution2.pvd')

# Time-stepping

t = 0

for n in range(num_steps):

# Update current time

t += dt

# Compute solution

solve(F == 0, u, bc)

# Save to file and plot solution

vtkfile << (u, t)

plot(u)

# Update previous solution

u_n.assign(u)

# Hold plot

interactive()

I have a term in my diffusion pde which represents the area where my function h is non-zero. I.E. ind(u) = 1 if u>=0 and ind(u) = 0 otherwise.

I have tried a couple of things and cannot figure out how to fit this into FEniCS. See the following:

def ind(u_n):

if u_n >=0:

return 1

else:

return 0

F = u*v*dx + dt*q(u)*dot(grad(u), grad(v))*dx - u_n*v*dx + ind(u_n)*(u_n)*v*dx

doesn't work, neither does

F = u*v*dx + dt*q(u)*dot(grad(u), grad(v))*dx - u_n*v*dx + E*u_n[u_n>0]*v*dx

Can anyone help me? Full code attached at the end.

Thanks, Zhen.

from fenics import *

import time

import scipy

def q(u):

"return nonlinear coefficient"

return u**3

T = 2.0 # final time

num_steps = 20 # number of time steps

dt = T / num_steps # time step size

E = 0.1 # evaporation constant

# Create mesh and define function space

nx = ny = 30

mesh = RectangleMesh(Point(-3, -3), Point(3, 3), nx, ny)

V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition

def boundary(x, on_boundary):

return on_boundary

bc = DirichletBC(V, Constant(0), boundary)

# Define initial value

u_0 = Expression('x[0]>=-1 && x[0]<=1 && x[1] >=-1 && x[1] <= 1 ? exp(-3*a*pow(x[0], 2) - a*pow(x[1], 2)):0',

degree=2, a=5)

#u_0 = Constant(0)

u_n = interpolate(u_0, V)

# Define variational problem

u = Function(V)

v = TestFunction(V)

f = Constant(0)

#note whether there is evaporation or not

def ind(u_n):

if u_n >=0:

return 1

else:

return 0

F = u*v*dx + dt*q(u)*dot(grad(u), grad(v))*dx - u_n*v*dx + E*u_n[u_n>0]*v*dx

# Create VTK file for saving solution

vtkfile = File('two_d_spreading/solution2.pvd')

# Time-stepping

t = 0

for n in range(num_steps):

# Update current time

t += dt

# Compute solution

solve(F == 0, u, bc)

# Save to file and plot solution

vtkfile << (u, t)

plot(u)

# Update previous solution

u_n.assign(u)

# Hold plot

interactive()

### 1 Answer

4

For your first approach consider this:

```
import numpy as np
def ind(u):
u_ind = Function(V)
pos = np.where(u.vector()>0.0)[0]
u_ind.vector()[pos] = 1.0
return u_ind.vector()
u_ind = Function(V)
u_ind.vector()[:] = ind(u_n)
```

And for the second one:

```
import numpy as np
u_ind = Function(V)
pos = np.where(u_n.vector()>0.0)[0]
u_ind.vector()[pos] = u_n.vector()[pos]
```

Thanks Hernan! May I ask how should I incorporate this into my F? Ultimately, I want to solve a pde where there is a term proportional to the area of the region where u is non zero.

best,

zhen

best,

zhen

written
10 months ago by
zhen shao

I was thinking in something like this:

```
u_ind = Function(V)
F = u*v*dx + dt*q(u)*dot(grad(u), grad(v))*dx - u_n*v*dx + E*u_ind*v*dx
# Create VTK file for saving solution
vtkfile = File('two_d_spreading/solution2.pvd')
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Compute solution
pos = np.where(u_n.vector()>0.0)[0]
u_ind.vector()[pos] = u_n.vector()[pos]
solve(F == 0, u, bc)
.
.
.
```

written
10 months ago by
Hernán Mella

That works! Thanks!

Thank you for helping me Mella, I am new to Python, do you recommend any online tutorial that I could go through?

Thank you for helping me Mella, I am new to Python, do you recommend any online tutorial that I could go through?

written
10 months ago by
zhen shao

Take a look here

Regards!

Regards!

written
10 months ago by
Hernán Mella

Please login to add an answer/comment or follow this question.