simple question about derivation

Question

Hi
I have a very simple question about taking derivative of an arbitrary function. For example :
g = x^2 + y^3

I want to take derivative of g w.r.t "x" and also "y" :

g="x[0]**x[0] + x[1]*x[1]*x[1]"
s = Expression(g)
derivation1 = s.dx(0)
derivation2 = s.dx(1)
print (derivation1)
print (derivation2)

But I am getting some errors. How do I have to write the code to get 2x after taking derivative w.r.t x or getting 3y^2 after taking derivative w.r.t y?

Answer 1

Expressions are interpolated into a finite element space before derivatives are evaluated. In particular, if you interpolate into a piecewise linear space all second derivatives will vanish. Considering the following example.

P1 = FiniteElement("Lagrange", triangle, 1)
P2 = FiniteElement("Lagrange", triangle, 2)

f = Expression("x[0]*x[1]", element = P1)
g = Expression("x[0]*x[1]", element = P2)


print assemble(f.dx(0).dx(1) * dx(domain = mesh))  # zero
print assemble(g.dx(0).dx(1) * dx(domain = mesh))  # nonzero

For simple functions of the coordinates, also consider the following alternative approach.

x, y = SpatialCoordinate(mesh)
h = x * y

print assemble(h.dx(0).dx(1) * dx)

This avoids interpolation into FEM space. Instead, the function is evaluated directly in the quadrature points.

Comments

Thank you for your response.
The thing you said gives back the integration of (d^2)h/dxdy= x*y over the domain. But what I am looking for is the value of derivative of a function at a specific point. For example:

from dolfin import *
import mshr
domain = mshr.Rectangle(Point(-1, -1), Point(2, 2))
mesh   = mshr.generate_mesh(domain, 20)
x,y = SpatialCoordinate(mesh)
h = x*y
a= h.dx(0)
print a(1,2)

In this case I want to get the value of the derivative of the h with respect to x (dh/dx = y) at point (1,2). I expect to get 2 but it gives back some error. Could you please help me fix it?
Thanks

---

You can project into a FunctionSpace before evaluating. For example,

V = FunctionSpace(P1, mesh)
a = project(h.dx(0), V)

print a(0.5, 0.5)

The projection should not be necessary, but it seems to be the easy way for now.

---

Thanks!
That was really helpful.