### Jump() in Expression over interface in DG Poisson

581

views

0

I have a question regarding how to implement a jump condition in FEniCS over an internal boundary. I want to use jump(u) within the expression for g (simplified in the code below). I want to use jump(u) on the interface marked '1' within an expression containing an exponential--in a manner similar to the expression for g below. Is this possible? If not, would it be better to extract the nodal values along interface '1', do my computation on the array? The error I get with g assigned as it is in the code below is "TypeError: expected default arguments for member variables to be scalars or GenericFunctions."

```
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
ofilename = 'model_axisymm2D.xml'
mesh = Mesh(ofilename)
subdomains = MeshFunction("size_t", mesh, "%s_physical_region.xml"%ofilename.split('.')[0])
boundaries = MeshFunction("size_t", mesh, "%s_facet_region.xml"%ofilename.split('.')[0])
V = FunctionSpace(mesh, 'DG', 2)
v = TestFunction(V)
u = TrialFunction(V)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dS = Measure('dS', domain=mesh, subdomain_data=boundaries)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
r = Expression('x[0]', degree=1, domain=mesh)
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2
alpha = Constant(10.0)
gamma = Constant(alpha)
VOLT = 40.0
g = Expression('alpha*exp(U*U)', alpha=1.0, U = jump(u), degree=2)
Fu_interior = inner(grad(v), grad(u))*r*dx
Fu_edges = -inner(jump(v,n), avg(grad(u)))*r*dS(0) \
- inner(avg(grad(v)), jump(u,n))*r*dS(0) \
+ alpha/h_avg*inner(jump(u), jump(v))*r*dS(0)
Fu_edges += - inner(avg(grad(v)), jump(u,n)-g*n('+'))*r*dS(1) \
+ alpha/h_avg*inner(jump(u)-g, jump(v))*r*dS(1)
Fu_bc = -inner(grad(v), (u-VOLT)*n)*r*ds(4) \
- inner(v*n, grad(u))*r*ds(4) \
+ (gamma/h)*v*(u-VOLT)*r*ds(4)
Fu_bc += -inner(grad(v), (u+VOLT)*n)*r*ds(6) \
- inner(v*n, grad(u))*r*ds(6) \
+ (gamma/h)*v*(u+VOLT)*r*ds(6)
Fu = Fu_interior + Fu_edges + Fu_bc
a = lhs(Fu)
l = rhs(Fu)
u = Function(V)
solve(a==l, u)
p = plot(u)
plt.colorbar(p)
```

Community: FEniCS Project

### 1 Answer

2

**Update:**I think I found the answer to my own problem here. Something to the tune of

`g = exp(1e-9*jump(u0)*jump(u0))`

works--I just have to use the solution from a previous timestep u0 to get my jump condition to work, which is not in the MWE I posted above.
1

I was playing around with a Boltzmann distribution, if it helps.

written
7 months ago by
Dan Sweeney

I see!

Thank you very much!

Thank you very much!

written
7 months ago by
Murilo Moreira

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

I was just curious about the physics of your problem.

What is it?

Thanks in advance.