Hello all
I am working on a tricky situation here. I need to use FEniCS’s automatic derivative to get gradient of a piecewise function. Suppose u(x,y)\in[0,1] is a solution of a 2D problem (I’m just using random_fun to make an example), f(u) is a 3-piece function based on the local values of u and f is actually continuously differentiable (see figure).
from dolfin import *
from numpy import array, polyval, random
import matplotlib.pyplot as plt
mesh = UnitSquareMesh(96, 96)
V = FunctionSpace(mesh, "CG", 1)
class random_fun(UserExpression):
def __init__(self, **kwargs):
super().__init__(**kwargs)
def eval(self, value, x):
value[0] = random.rand()
def value_shape(self):
return ()
u = interpolate( random_fun(), V)
plt.figure()
fig=plot(u)
plt.colorbar(fig)
plt.show()
class f(UserExpression):
def __init__(self, u_fun, **kwargs):
super().__init__(**kwargs)
self.u_fun = u_fun
def eval(self, value, x):
if self.u_fun(x) <= 0.25/2.:
value[0] = self.u_fun(x)/0.25
elif self.u_fun(x) >= 2.*0.25:
value[0] = 1.
else:
value[0] = polyval(array([256./27., -128./9., 64./9., -5./27.]), self.u_fun(x))
def value_shape(self):
return ()
fu = project(f(u),V)
plt.figure()
fig=plot(fu)
plt.colorbar(fig)
plt.show()
du = TestFunction(V)
f_form = fu * dx
dfdu = derivative(f_form, u, du)
dfdu_fun = Function(V, assemble(dfdu))
plt.figure()
fig=plot(dfdu_fun)
plt.colorbar(fig)
plt.show()
I also tried to plot u, f(u), and \frac{df(u)}{du} to visualize and verify but wasn’t successful. I am using FEniCS version 2019.1.
My understanding of your question is as follows: you have a function F(x, y) = f(u(x, y)), where f(w) is a scalar function of a scalar input w, and f is defined as shown in the figure in your post. When you say you want to get the gradient of the function, do you mean that you want:
- The gradient, i.e. \nabla F = \frac{\partial F}{\partial x} \hat{i} + \frac{\partial F}{\partial y} \hat{j}
- The derivative F'(x, y) = f'(u(x, y)), where f'(w) = \frac{df(w)}{dw}
- The functional derivative (a.k.a. the Gateaux derivative) of the functional G, which it appears you are trying to calculate in your code (because of the use of
derivative):
dG(u; \delta u) = \lim_{\tau \rightarrow 0} {\frac{d}{d \tau} G(u + \tau \, \delta u)} \\
G(u) = \int_{[0, 1] \times [0, 1]} {f(u(x, y)) \, \mathrm{d}x}
2 Likes
What do you think about this?
from dolfin import *
from numpy import array, polyval, random
import matplotlib.pyplot as plt
mesh = UnitSquareMesh(96, 96)
V = FunctionSpace(mesh, "CG", 1)
class random_fun(UserExpression):
def __init__(self, **kwargs):
super().__init__(**kwargs)
def eval(self, value, x):
value[0] = random.rand()
def value_shape(self):
return ()
u = interpolate( random_fun(), V)
plt.figure()
fig=plot(u)
plt.colorbar(fig)
plt.show()
fu = conditional( u<= 0.25/2., u/0.25, conditional( u >= 2.*0.25, 1., 256./27.*u**3. -128./9.*u**2. + 64./9.*u -5./27. ) )
plt.figure()
fig=plot(fu)
plt.colorbar(fig)
plt.show()
du = TestFunction(V)
f_form = fu * dx
dfdu = derivative(f_form, u, du)
dfdu_fun = Function(V, assemble(dfdu))
plt.figure()
fig=plot(dfdu_fun)
plt.colorbar(fig)
plt.show()
That appears to be correct. To test for correctness, I’d suggest using a small mesh (e.g. mesh = UnitSquareMesh(1, 1)), assembling (e.g. b = assemble(dfdu)) and checking the entries of the vector to see if they are correct for simple cases. E.g.
u.vector()[:] = 0.1 # df/du = 4 for all x
b = assemble(dfdu)
print(b[:]) # should give [2/3, 4/3, 4/3, 2/3]
u.vector()[:] = 0.75 # df/du = 0 for all x
b = assemble(dfdu)
print(b[:]) # should give [0., 0., 0., 0.]
u.vector()[:] = 0.25 # calculate m = df/du
b = assemble(dfdu)
print(b[:]) # should give [m / 6, m / 3, m / 3, m / 6]
2 Likes
