# Is there a way to get the partial derivative for a gradient vector?

**URL:** <https://fenicsproject.discourse.group/t/is-there-a-way-to-get-the-partial-derivative-for-a-gradient-vector/15753>\
**Category:** Legacy DOLFIN\
**Tags:** mesh, ufl\
**Created:** [September 9, 2024, 11:55am UTC](https://fenicsproject.discourse.group/t/is-there-a-way-to-get-the-partial-derivative-for-a-gradient-vector/15753 "2024-09-09T11:55:39Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Swani](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swani/32/3940_2.png) [@Swani](https://fenicsproject.discourse.group/u/Swani)\
**Post date:** [September 9, 2024, 11:55am UTC](https://fenicsproject.discourse.group/t/is-there-a-way-to-get-the-partial-derivative-for-a-gradient-vector/15753/1 "2024-09-09T11:55:39Z")

</div>

Hello. Dear community 👋

I was wondering if it’s possible **to get the partial derivative(Cartesian coordinate) for a gradient vector in legacy FEniCS**.  
(this might be a bit silly, but I’m posting it for the sake of curiosity, please forgive me)

I was trying to implement the following formula in legacy FEniCS.

u\_{elem}(x,y)=x^2+y^2 \to \nabla{u\_{elem}(x,y)}=\left[{2x\atop2y}\right] \to {\partial{\nabla{u\_{elem}(x,y)}}\over\partial{x}}=\left[{2\atop0}\right], {\partial{\nabla{u\_{elem}(x,y)}}\over\partial{y}}=\left[{0\atop2}\right]

The code looks like:

```auto
from dolfin import *

dom = RectangleMesh(Point(0.0, 0.0), Point(2.0, 2.0), 2, 2, diagonal="right")
Vh = VectorFunctionSpace(dom, "CG", 1)
CG1 = FunctionSpace(dom, "CG", 1)
u_nodal = interpolate(Expression("x[0]*x[0]+x[1]*x[1]", degree=2), CG1)
print(f"size of u_nodal vector: {u_nodal.vector().get_local().size}")
print(f"u_nodal vector: {u_nodal.vector().get_local()}\n")

grd_u_nodal = project(grad(u_nodal), Vh)
print(f"size of grad(u_nodal) vector: {grd_u_nodal.vector().get_local().size}")
print(f"grad(u_nodal) vector: {grd_u_nodal.vector().get_local()}\n")

x0, x1 = MeshCoordinates(dom)
print(f"dgrad(u_nodal)dx: {assemble(diff(grad(u_nodal), variable(x0))*dx)}")

```

The output and error:

```auto
size of u_nodal vector: 9
u_nodal vector: [4. 1. 5. 0. 2. 8. 1. 5. 4.]

size of grad(u_nodal) vector: 18
grad(u_nodal) vector: [0.85714286 3.14285714 0.57142857 2.28571429 1.71428571 3.42857143
 0.28571429 0.28571429 2. 2. 3.71428571 3.71428571
 2.28571429 0.57142857 3.42857143 1.71428571 3.14285714 0.85714286]

```

```auto
Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ().

```

I’m thinking that the error is probably caused by a difference between the array shapes, **is there any way to fix this?**  
**(And the reason for the decimals in the gradient array is numerical difference?)**

I hope this clears up some of my questions, however minor. 🙏

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [September 9, 2024, 4:05pm UTC](https://fenicsproject.discourse.group/t/is-there-a-way-to-get-the-partial-derivative-for-a-gradient-vector/15753/2 "2024-09-09T16:05:21Z")

</div>

Anything wrong with `u_nodal.dx(0).dx(0)`?

Of course, you’d better be using a higher degree than 1 in your `Vh`, otherwise the second derivative will surely be zero.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 9, 2024, 7:11pm UTC](https://fenicsproject.discourse.group/t/is-there-a-way-to-get-the-partial-derivative-for-a-gradient-vector/15753/3 "2024-09-09T19:11:52Z")

</div>

You cannot assemble `vector`-expressions by themselfs. In Legacy Dolfin (and DOLFINx).  
You can print `grd_u_nodal` because you have used project, that creates an inner product for the `VectorFunctionSpace`.

Please note that you can use

```python
import dolfin

mesh = dolfin.UnitSquareMesh(10, 10)

x, y = dolfin.SpatialCoordinate(mesh)

u = x **2 + y** 2

grad_u = dolfin.grad(u)

grad_u_ex = dolfin.as_vector((2 * x, 2 * y))

error = dolfin.assemble(
    dolfin.inner(grad_u - grad_u_ex, grad_u - grad_u_ex) * dolfin.dx
)
print(error)

Hessian_u = dolfin.grad(grad_u)

Hessian_u_ex = 2 * dolfin.Identity(2)

error = dolfin.assemble(
    dolfin.inner(Hessian_u - Hessian_u_ex, Hessian_u - Hessian_u_ex) * dolfin.dx
)
print(error)

```

which gives you exact accuracy, for instance if u\_elem = x **3 + y** 3 then

```python
import dolfin

mesh = dolfin.UnitSquareMesh(10, 10)

x, y = dolfin.SpatialCoordinate(mesh)

u = x **3 + y** 3

grad_u = dolfin.grad(u)

grad_u_ex = dolfin.as_vector((3 * x **2, 3 * y** 2))

error = dolfin.assemble(
    dolfin.inner(grad_u - grad_u_ex, grad_u - grad_u_ex) * dolfin.dx
)
print(error)

Hessian_u = dolfin.grad(grad_u)

Hessian_u_ex = dolfin.as_tensor(((6 * x, 0), (0, 6 * y)))
error = dolfin.assemble(
    dolfin.inner(Hessian_u - Hessian_u_ex, Hessian_u - Hessian_u_ex) * dolfin.dx
)
print(error)

```

gives you

```bash
0.0
0.0

```

---

<div class="post-metadata">

**Author:** ![Swani](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swani/32/3940_2.png) [@Swani](https://fenicsproject.discourse.group/u/Swani)\
**Post date:** [September 10, 2024, 4:43am UTC](https://fenicsproject.discourse.group/t/is-there-a-way-to-get-the-partial-derivative-for-a-gradient-vector/15753/4 "2024-09-10T04:43:57Z")

</div>

Many responses! 🙌  
Thank you for always being available to reply my post. You’re a great help in solving my curiosity!

Dear @francesco-ballarin,  
From comments, what does this mean?

> [@francesco-ballarin](#):
>
> u\_nodal.dx(0).dx(0)

The output looked like that.  
(I assume it was outputting a UFL form, am I right?)

```auto
(grad((grad(f_1988))[0]))[0]

```

Dear @dokken,  
In reply, you used `grad(u)` and `grad_u_ex` to show that there was no difference between the two parameters. It was an informative(as was the Hessian code result also).

Furthermore, while checking your reply, I have some additional questions.  
a. In the code, you showed the Hessian, shouldn’t the Hessian( \nabla^{2}u\_{elem}(x,y) ) be viewed differently from the partial derivative of the gradient( {\partial\nabla{u\_{elem}(x,y)}\over\partial{x} \ or \ \partial{y}} )?  
(If my knowledge is wrong, please correct me.)

b. Since `grad(u)` is a UFL form, vector array extraction seems to be difficult, so I tried vector array extraction with `as_vector`, but with difficulty. Is there a way to extract either `grad_u` or `grad_u_ex` as a vector array? As far as I know, there doesn’t seem to be a way.  
While this is an important example for noticing that the two results are not different, I think there is a difficulty in utilising the gradient(or Hessian) array.  
(The reason I want to extract the values into a vector array is to compare them later with other operations).

c. This is a minor question, but is there a specific reason why you used inner product norm(similar as L^{2} norm squared) in the error checking?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 10, 2024, 5:39am UTC](https://fenicsproject.discourse.group/t/is-there-a-way-to-get-the-partial-derivative-for-a-gradient-vector/15753/5 "2024-09-10T05:39:56Z")

</div>

> [@Swani](#):
>
> Furthermore, while checking your reply, I have some additional questions.  
> a. In the code, you showed the Hessian, shouldn’t the Hessian( \nabla^{2}u\_{elem}(x,y) ∇2uelem(x,y) \nabla^{2}u\_{elem}(x,y) ) be viewed differently from the partial derivative of the gradient( {\partial\nabla{u\_{elem}(x,y)}\over\partial{x} \ or \ \partial{y}} ∂∇uelem(x,y)∂x or ∂y{\partial\nabla{u\_{elem}(x,y)}\over\partial{x} \ or \ \partial{y}} )?  
> (If my knowledge is wrong, please correct me.)

You can do that by calling`grad_u.dx(0)` or `grad_u.dx(1)`.

> [@Swani](#):
>
> b. Since `grad(u)` is a UFL form, vector array extraction seems to be difficult, so I tried vector array extraction with `as_vector`, but with difficulty. Is there a way to extract either `grad_u` or `grad_u_ex` as a vector array? As far as I know, there doesn’t seem to be a way.  
> While this is an important example for noticing that the two results are not different, I think there is a difficulty in utilising the gradient(or Hessian) array.  
> (The reason I want to extract the values into a vector array is to compare them later with other operations).

Please be clear on how you want to do this comparison.  
As the examples you have shown has been simple, analytic expressions, you could simply interpolate the analytical expression into the function space you want to work with. See for instance

> [@Defining real part of complex equation using Expression() in FEniCS](https://fenicsproject.discourse.group/t/defining-real-part-of-complex-equation-using-expression-in-fenics/15631/2):
>
> Legacy dolfin doesn’t really have complex support. You can use sympy for this though: from dolfin import \* mesh = UnitSquareMesh(10, 10) import sympy x, y, t = sympy.symbols("x[0], x[1] t", real=True) v\_in = y \* (1 - y) \* sympy.exp(1j \* t) v\_real, \_ = v\_in.as\_real\_imag() expr = Expression(str(v\_real), degree=2, t=0.0) print(assemble(expr \* dx(domain=mesh))) expr.t = 2 print(assemble(expr \* dx(domain=mesh)))

> [@Convert sympy to ufl?](https://fenicsproject.discourse.group/t/convert-sympy-to-ufl/3852/2):
>
> Using strings and eval: import sympy x = sympy.Symbol("x") cos = sympy.cos sympy\_input = cos(x) + x \*\* 2 import fenics mesh = fenics.UnitIntervalMesh(10) cos = fenics.cos x = fenics.SpatialCoordinate(mesh)[0] ufl\_output = eval(str(sympy\_input)) print(fenics.assemble(ufl\_output\*fenics.dx)) However, is there a specfic reason for using sympy? UFL-forms can be differentiated etc, so in most use cases there is no need for using sympy. EDIT: You can also use: [https://github.com/MiroK/ulfy](https://github.com/MiroK/ulfy)

> [@Swani](#):
>
> c. This is a minor question, but is there a specific reason why you used inner product norm(similar as L^{2} L2L^{2} norm squared) in the error checking?

As pointed out in the previous part of the question, what I have described gives you a symbolic expression that you can evaluate when integrating. Thus I chose an L2 norm, as it Seems natural in this setting

---

<div class="post-metadata">

**Author:** ![Swani](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swani/32/3940_2.png) [@Swani](https://fenicsproject.discourse.group/u/Swani)\
**Post date:** [September 11, 2024, 5:59am UTC](https://fenicsproject.discourse.group/t/is-there-a-way-to-get-the-partial-derivative-for-a-gradient-vector/15753/6 "2024-09-11T05:59:40Z")

</div>

Thank you for comment, @dokken! It was very helpful. I really appreciate it!  
_(Apologies for the late reply. 🙏)_

In case it helps anyone reading this in the future, here’s my understanding based on the code you answered(bottom of my reply, “MWE”).  
To wrap things up, I noticed a difference in the `u` array results between the code I originally attached and the code you replied to, and between **projection** and **interpolation**.

Despite being a simple mesh(I changed the mesh size 10 \times 10 to 2 \times 2 to make it easier to see).

Code:

```auto
from dolfin import *

mesh = UnitSquareMesh(2, 2)
CG1 = FunctionSpace(mesh, "CG", 1)

x, y = SpatialCoordinate(mesh)
u = x **2 + y** 2
print(f"size of u vector: {project(u, CG1).vector().get_local().size}")
print(f"u vector: {project(u, CG1).vector().get_local()}")
print(f"interpolate u: {interpolate(Expression("x[0]*x[0]+x[1]*x[1]", degree=2), CG1).vector().get_local()}\n")

```

Output:

```auto
size of u vector: 9
u vector: [ 0.93571429 0.16428571 1.16428571 -0.09285714 0.42142857 1.90714286
  0.16428571 1.16428571 0.93571429]
interpolate u: [1. 0.25 1.25 0. 0.5 2. 0.25 1.25 1.]

```

There was a difference.  
(When I say “gradient operation” later, I mean the same thing as in the code I referenced, but only interpolating the `u` part, or running it with a difference of `u=x **2+y** 2`)

For me, intuitively, the interpolated result of the two values looks right, but when I did a “gradient operation” with `interpolate u`, it was different from the exact solution.  
On the other hand, when I checked the `u` array through `project`, it intuitively looked like there was an error, but when I did the “gradient operation”, it matched the exact solution.

In my opinion, defining `u` as `u=x **2+y** 2` rather than using `interpolate` seems to be wrong to the eye, but I think it is correct considering the operation, **is there any way to improve this**?  
_(I was hoping to be able to use projection and still see results that were reasonably consistent with the interpolation results)_

MWE:

```auto
from dolfin import *

mesh = UnitSquareMesh(2, 2)
CG1 = FunctionSpace(mesh, "CG", 1)
Vh = VectorFunctionSpace(mesh, "CG", 1)
Th = TensorFunctionSpace(mesh, "CG", 1)

x, y = SpatialCoordinate(mesh)
u = x **2 + y** 2

# Hmm... different!
print(f"size of u vector: {project(u, CG1).vector().get_local().size}")
print(f"u vector: {project(u, CG1).vector().get_local()}")
print(f"interpolate u: {interpolate(Expression("x[0]*x[0]+x[1]*x[1]", degree=2), CG1).vector().get_local()}\n")

grad_u = grad(u)
grad_u_ex = as_vector((2*x, 2*y))

print(f"grad_u size: {project(grad_u, Vh).vector().size()}")
print(f"grad_u vector: {project(grad_u, Vh).vector()[:]}")
print(f"x of grad_u size: {project(grad_u.dx(0), Vh).vector().get_local().size}")
print(f"x of grad_u vector: {project(grad_u.dx(0), Vh).vector().get_local()}")
print(f"y of grad_u size: {project(grad_u.dx(1), Vh).vector().get_local().size}")
print(f"y of grad_u vector: {project(grad_u.dx(1), Vh).vector().get_local()}")
grad_error = assemble(inner(grad_u - grad_u_ex, grad_u - grad_u_ex) * dx)
print(f"error between grad_u and grad_u_ex: {grad_error}\n")

Hessian_u = grad(grad_u)
Hessian_u_ex = 2*Identity(2)

print(f"Hessian size: {project(Hessian_u, Th).vector().get_local().size}")
print(f"Hessian vector: {project(Hessian_u, Th).vector().get_local()}")
hes_error = assemble(inner(Hessian_u - Hessian_u_ex, Hessian_u - Hessian_u_ex) * dx)
print(f"error between Hessian_u and Hessian_u_ex: {hes_error}")

```

Output:

```auto
size of u vector: 9
u vector: [ 0.93571429 0.16428571 1.16428571 -0.09285714 0.42142857 1.90714286
  0.16428571 1.16428571 0.93571429]
interpolate u: [1. 0.25 1.25 0. 0.5 2. 0.25 1.25 1.]

grad_u size: 18
grad_u vector: [-9.58434720e-17 2.00000000e+00 2.51534904e-17 1.00000000e+00
  1.00000000e+00 2.00000000e+00 8.32667268e-17 6.44227717e-17
  1.00000000e+00 1.00000000e+00 2.00000000e+00 2.00000000e+00
  1.00000000e+00 7.13714802e-17 2.00000000e+00 1.00000000e+00
  2.00000000e+00 -8.32667268e-17]
x of grad_u size: 18
x of grad_u vector: [2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0.]
y of grad_u size: 18
y of grad_u vector: [0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2.]
error between grad_u and grad_u_ex: 0.0

Hessian size: 36
Hessian vector: [2. 0. 0. 2. 2. 0. 0. 2. 2. 0. 0. 2. 2. 0. 0. 2. 2. 0. 0. 2. 2. 0. 0. 2.
 2. 0. 0. 2. 2. 0. 0. 2. 2. 0. 0. 2.]
error between Hessian_u and Hessian_u_ex: 0.0

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 11, 2024, 7:26am UTC](https://fenicsproject.discourse.group/t/is-there-a-way-to-get-the-partial-derivative-for-a-gradient-vector/15753/7 "2024-09-11T07:26:10Z")

</div>

> [@Swani](#):
>
> To wrap things up, I noticed a difference in the `u` array results between the code I originally attached and the code you replied to, and between **projection** and **interpolation**.
> 
> Despite being a simple mesh(I changed the mesh size 10 \times 10 10×1010 \times 10 to 2 \times 2 2×22 \times 2 to make it easier to see).
> 
> Code:
> 
> ```auto
> from dolfin import *
> 
> mesh = UnitSquareMesh(2, 2)
> CG1 = FunctionSpace(mesh, "CG", 1)
> 
> x, y = SpatialCoordinate(mesh)
> u = x **2 + y** 2
> print(f"size of u vector: {project(u, CG1).vector().get_local().size}")
> print(f"u vector: {project(u, CG1).vector().get_local()}")
> print(f"interpolate u: {interpolate(Expression("x[0]*x[0]+x[1]*x[1]", degree=2), CG1).vector().get_local()}\n")
> 
> ```
> 
> Output:
> 
> ```auto
> size of u vector: 9
> u vector: [ 0.93571429 0.16428571 1.16428571 -0.09285714 0.42142857 1.90714286
> 0.16428571 1.16428571 0.93571429]
> interpolate u: [1. 0.25 1.25 0. 0.5 2. 0.25 1.25 1.]
> 
> ```
> 
> There was a difference.

Of course there will be a difference between interpolation and projection, see for instance: [Lecture 8 — Lecture notes MATMEK-4270](https://matmek-4270.github.io/matmek4270-book/lecture8.html)

> [@Swani](#):
>
> In my opinion, defining `u` as `u=x **2+y** 2` rather than using `interpolate` seems to be wrong to the eye, but I think it is correct considering the operation, **is there any way to improve this**?

Using `Expression` the string describing `u` will be evaluated at interpolation points, while using `x, y = SpatialCoordinate(mesh)` gives you an evaluation at the quadrature points used in the projection scheme.  
This means that you can control the accuracy of approximation by choosing a suitable quadrature rule (up to the best possible approximation).

---

<div class="post-metadata">

**Author:** ![Swani](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swani/32/3940_2.png) [@Swani](https://fenicsproject.discourse.group/u/Swani)\
**Post date:** [September 12, 2024, 3:16am UTC](https://fenicsproject.discourse.group/t/is-there-a-way-to-get-the-partial-derivative-for-a-gradient-vector/15753/8 "2024-09-12T03:16:32Z")

</div>

Thanks for your reply!  
I’ll do some more studying based on what you’ve show me! Thanks! 😁  
I’ll leave a comment if I don’t understand something later!