$\nabla.(f(x)\vec{r})$ to interpolate, or not to interpolate

Question

I am trying to solve variable coefficient Poisson equation in mixed formulation using periodic boundary conditions.

$\sigma-\nabla(\phi)=0$
$\nabla.(\epsilon(r)\sigma)=4\pi\rho(r)$

until today, I was able to get some results using the following approach

Eps2=1.0
Eps1=2.45
Z0_B=1.78
Z0_T=6.22
Beta=0.208645
eps_str="0.5*({eps2}-{eps1})*erf((x[2]-{z0_B})/{beta})*erf((x[2]-{z0_T})/{beta}) + 0.5*({eps1}+{eps2})"\
        .format(eps2=Eps2,eps1=Eps1,z0_B=Z0_B,z0_T=Z0_T,beta=Beta)
epsilon=Expression(eps_str)
a = (dot(sigma, tau) + div(tau)*u + div(epsilon*sigma)*v)*dx - inner(u, nu)*dx - inner(mu, v)*dx

However, when I was playing with some other $\epsilon$ models, something broke, and I started getting

Cannot get geometric dimension from an expression with no domains!

I am puzzled why it was working before, and why it stopped working now. My immediate fix was to replace epsilon to an interpolation to a DG function space

DG = FunctionSpace(mesh, "DG", 0, constrained_domain=pb)
epsilon=interpolate(Expression(eps_str),DG)

but I am not so sure if this is a mistake or not. Especially since I am using an adaptive mesh algorithm.

Comments

Another attempt by taking the gradient by hand results in a NAN

eps_str="0.5*({eps2}-{eps1})*erf((x[2]-{z0_B})/{beta})*erf((x[2]-{z0_T})/{beta}) + 0.5*({eps1}+{eps2})"\
           .format(eps2=Eps2,eps1=Eps1,z0_B=Z0_B,z0_T=Z0_T,beta=Beta)
eps_gradx_str="0"
eps_grady_str="0"
eps_gradz_str="0.56418958355*(({eps2}-{eps1})/{beta})*(exp(-pow((x[2]-{z0_B})/{beta},2))*erf((x[2]-{z0_T})/{beta})+erf((x[2]-{z0_B})/{beta})*exp(-pow((x[2]-{z0_T})/{beta},2)))"
          .format(eps2=Eps2,eps1=Eps1,z0_B=Z0_B,z0_T=Z0_T,beta=Beta)
BDM = FunctionSpace(mesh, "BDM", 1, constrained_domain=pb)
DG = FunctionSpace(mesh, "DG", 0, constrained_domain=pb)
R = FunctionSpace(mesh, 'R', 0, constrained_domain=pb)
W = MixedFunctionSpace([BDM, DG, R])
(sigma, u, mu) = TrialFunctions(W)
(tau, v, nu) = TestFunctions(W)
epsilon=Expression(eps_str,element=DG.ufl_element())
epsilon_grad=Expression((eps_gradx_str,eps_grady_str,eps_gradz_str),element=BDM.ufl_element()) #now it is a vector
a = (dot(sigma, tau) + div(tau)*u + epsilon*v*div(sigma)+inner(epsilon_grad,sigma)*v)*dx - inner(u, nu)*dx - inner(mu, v)*dx

When I look at how the epsilon_grad looks by interpolating via a DGM function, the vectors are not in z and -z directions as expected. X and Y components are also present. Can this be the problem?

Answer 1

Try:

epsilon = Expression(eps_str, domain=mesh)

and I think you should be fine.

Comments

epsilon=Expression(eps_str,element=DG.ufl_element())
seems to be working as well.