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

581
views
0
7 months ago by

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)

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

+ (gamma/h)*v*(u-VOLT)*r*ds(4)

+ (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

2
7 months ago by
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.
Dear Dan,