# How to Implement First and Second Derivatives w.r.t. Specific Directions in FEniCS
**URL:** https://fenicsproject.discourse.group/t/how-to-implement-first-and-second-derivatives-w-r-t-specific-directions-in-fenics/16536
**Category:** General
**Tags:** dolfinx
**Created:** [December 17, 2024, 4:49pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-first-and-second-derivatives-w-r-t-specific-directions-in-fenics/16536 "2024-12-17T16:49:05Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@Great\_Mammuth](https://fenicsproject.discourse.group/u/Great_Mammuth)
#### Post date: [December 17, 2024, 4:49pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-first-and-second-derivatives-w-r-t-specific-directions-in-fenics/16536/1 "2024-12-17T16:49:05Z")
Dear everyone,
I am currently working on solving a PDE with first and second derivatives explicitly with respect to two spatial variables π₯1 and π₯2. I am having some difficulty identifying the correct FEniCS operators to implement this equation properly. The equation I am trying to solve is:
(πΉ,22 β β1) π,11 + (πΉ,11 β β2) π,22 β 2 πΉ,12 π,12 β π,1 π1 πβ π,2 π2 π β π = 0
Where
π(x1, x2) is the unknown function
π,i represents the first derivatives with respect to π₯π
π,ij represents the second derivatives with respect to π₯π and π₯π.
F(x1, x2) represents a scalar function of spatial variables x1 and x2 with its second variables.
h1, h2, π1 and π2 are scalars and p is a known load term.
I attempted to implement the equation using grad(f) and inner() terms, but this approach includes derivatives in all spatial directions, which does not match the PDE that requires derivatives explicitly with respect to π₯1 and π₯2 only.
My weak form currently looks like this:
```auto
a = (F_S22_minus_h1 * inner(grad(f), grad(v)) * dx
+ F_S11_minus_h2 * inner(grad(f), grad(v)) * dx
- 2 * F_S12 * inner(grad(f), grad(v)) * dx)
-load_times_lambda_1 * v * dx - + load_times_lambda_2 * v * dx
L = -load * v * dx
```
How can I explicitly apply first and second derivatives with respect to π₯1 and π₯2 only in the weak form of my PDE? Is there a better operator or method in FEniCS to achieve this? I apologize if this is a basic questionβI am still quite new to FEniCS and learning its syntax. Any guidance, hints, or references would be greatly appreciated!
Thank you very much, and I look forward to hearing from you.
Best Regards,
Luigi
---
### Author:  [@Stein](https://fenicsproject.discourse.group/u/Stein)
#### Post date: [December 17, 2024, 6:39pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-first-and-second-derivatives-w-r-t-specific-directions-in-fenics/16536/2 "2024-12-17T18:39:45Z")
Fenics permits simply `f.dx(0)` to take the derivative of a `Function` object `f` wrt x0. You can repeat this operation, ie. `f.dx(0).dx(0)`.
Note that for C0 continuous functions this is a βdangerousβ operation, and often is not what you actually wantβ¦ (Taking second derivatives, that is)
---
### Author:  [@Great\_Mammuth](https://fenicsproject.discourse.group/u/Great_Mammuth)
#### Post date: [December 23, 2024, 10:14pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-first-and-second-derivatives-w-r-t-specific-directions-in-fenics/16536/3 "2024-12-23T22:14:21Z")
Hi Stein,
Thank you very much for your advice. I have implemented your suggestions and compared the result with the one obtained by solving the same equation in Mathematica. In particular, I have simplified the original equation by setting lambda and h equal to zero as in the following code:
```auto
def PDE_solver_fenics(case, mesh_file, z_values, params, order, q):
# Load the mesh and define the function space
mesh = Mesh()
with XDMFFile(mesh_file) as infile:
infile.read(mesh)
V = FunctionSpace(mesh, "P", 2)
# Define trial and test functions
f = TrialFunction(V)
v = TestFunction(V)
# Define parameters and symbolic expressions
l, b = 2 * np.max(mesh.coordinates()[:, 0]), 2 * np.max(mesh.coordinates()[:, 1])
x1, x2 = sp.symbols('x1 x2')
A_expr = Func(case, b, l, x1, x2, params, order)
F_S11_expr = sp.diff(A_expr, x1, x1)
F_S22_expr = sp.diff(A_expr, x2, x2)
F_S12_expr = sp.diff(A_expr, x1, x2)
# Convert to numerical functions
F_S11_func = sp.lambdify((x1, x2), F_S11_expr, 'numpy')
F_S22_func = sp.lambdify((x1, x2), F_S22_expr, 'numpy')
F_S12_func = sp.lambdify((x1, x2), F_S12_expr, 'numpy')
# Define function terms
F_S11 = FuncTerm(F_S11_func, degree=2)
F_S22 = FuncTerm(F_S22_func, degree=2)
F_S12 = FuncTerm(F_S12_func, degree=2)
load = Constant(q)
a = (F_S22 * inner(grad(f), grad(v)) * dx
+ F_S11 * inner(grad(f), grad(v)) * dx
- 2 * F_S11 * inner(grad(f), grad(v)) * dx)
# a = (F_S22 * f.dx(0).dx(0) * v * dx +
# F_S11 * f.dx(1).dx(1) * v * dx -
# 2 * F_S12 * f.dx(0).dx(1) * v * dx)
# Linear form
L = -load * v * dx
# Apply boundary conditions and solve
bc = DirichletBC(V, Constant(0.0), "on_boundary")
f_solution = Function(V)
solve(a == L, f_solution, bc)
return f_solution, mesh
```
When using the terms inner and grad, the solution is similar to the one obtained from Mathematica. However, when I use the terms f.dx(0).dx(0) (commented in the code), the solution is not stable. Could this be due to the direct computation of second derivatives, as you mentioned earlier? Am I implementing these derivatives correctly in FEniCS?
Thanks again for your support, and I look forward to any further guidance you can provide! Also, Merry Christmas!
Best Regards,
Luigi
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [December 23, 2024, 10:30pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-first-and-second-derivatives-w-r-t-specific-directions-in-fenics/16536/4 "2024-12-23T22:30:16Z")
> [@Great\_Mammuth](#):
>
> When using the terms inner and grad, the solution is similar to the one obtained from Mathematica. However, when I use the terms f.dx(0).dx(0)
One doesnβt use the strong form of a PDE in the finite element method. One of the key aspects of FEM is that you set up a weak formulation (variational form), where you lessen the requirement of the function spaces used for your unknown). This has for instance been discussed in
> [@Is it possible to solve original or weak formulations of the boundary value problems in DOLFINx?](https://fenicsproject.discourse.group/t/is-it-possible-to-solve-original-or-weak-formulations-of-the-boundary-value-problems-in-dolfinx/16282/4):
>
> Biharmonic equation: [Biharmonic equation β DOLFINx 0.10.0.0 documentation](https://docs.fenicsproject.org/dolfinx/main/python/demos/demo_biharmonic.html) Honestly, I canβt give you a full review of Galerkin based finite element methods. I think a good source of information is: [pde - What is the purpose of using integration by parts in deriving a weak form for FEM discretization? - Computational Science Stack Exchange](https://scicomp.stackexchange.com/questions/7845/what-is-the-purpose-of-using-integration-by-parts-in-deriving-a-weak-form-for-fe#:~:text=bearing%20the%20same%20spirit%20with,not%20require%20integration%20by%20parts). Yes, if you want higher order regularity, you would need to use hermitian elements or splines. You say that it is Β«as easy asΒ». Sure, if someone has timβ¦
referring to [pde - What is the purpose of using integration by parts in deriving a weak form for FEM discretization? - Computational Science Stack Exchange](https://scicomp.stackexchange.com/questions/7845/what-is-the-purpose-of-using-integration-by-parts-in-deriving-a-weak-form-for-fe#:~:text=bearing%20the%20same%20spirit%20with,not%20require%20integration%20by%20parts)
---
### Author:  [@Great\_Mammuth](https://fenicsproject.discourse.group/u/Great_Mammuth)
#### Post date: [December 27, 2024, 5:45pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-first-and-second-derivatives-w-r-t-specific-directions-in-fenics/16536/5 "2024-12-27T17:45:43Z")
Hi dokken,
Thank you very much for your valuable comment. I have read the links you provided, however I was not able to understand how to properly write the PDE I want to solve in the weak form. Can you please have a quick look at the following equation and tell me if I have written it correctly?
The equation I am trying to solve is (itβs the same as reported at the beginning of the post for lambda=0):
πΉ,22 \* π,11 + πΉ,11 \* π,22 β 2 \* πΉ,12 \* π,12 β π = 0
Where
π(x1, x2): the unknown function to be solved for.
π,ij: the second partial derivatives with respect to xi and xj.
F(x1, x2) represents a scalar function of spatial variables x1 and x2.
F,ij the second partial derivatives with respect to xi and xj.
I have written the weak form of the equation as:
```auto
a = (F_S22 * f.dx(0).dx(0) * v * dx
+ F_S11 * f.dx(1).dx(1) * v * dx
- 2 * F_S12 * f.dx(0).dx(1) * v * dx)
```
Are the terms `v * dx` correct?
Thank you very much in advance!
Best Regards,
Luigi
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [December 27, 2024, 8:08pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-first-and-second-derivatives-w-r-t-specific-directions-in-fenics/16536/6 "2024-12-27T20:08:55Z")
They look fine to me, but again, one avoids second derivatives in standard Galerkin problems, thus you should integrate these terms by parts.
---
### Author:  [@Great\_Mammuth](https://fenicsproject.discourse.group/u/Great_Mammuth)
#### Post date: [December 28, 2024, 7:22pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-first-and-second-derivatives-w-r-t-specific-directions-in-fenics/16536/7 "2024-12-28T19:22:08Z")
Many thanks for your help so far! I have integrated by parts to eliminate the second derivatives in my equation. However, I am unsure how to handle the boundary terms that appear after integration by parts.
In particular, when integrating one of the terms:
\int\_\Omega F\_{,22} f\_{,11} v \, dx = -\int\_\Omega F\_{,22} f\_{,1} v\_{,1} \, dx + \big[F\_{,22} f\_{,1} v \big]\_{\partial \Omega}.
I have tried to omit the term \big[F\_{,22} f\_{,1} v \big]\_{\partial \Omega} and the solutions seems stable. However, I am unsure if this is the correct approach. Should I include this term explicitly? If so, how would I incorporate it correctly into the weak form?
Thank you very much for your assistance, it is greatly appreciated!
Regards,
Luigi
---
### Author:  [@Great\_Mammuth](https://fenicsproject.discourse.group/u/Great_Mammuth)
#### Post date: [December 29, 2024, 11:11pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-first-and-second-derivatives-w-r-t-specific-directions-in-fenics/16536/8 "2024-12-29T23:11:12Z")
A brief follow up on my previous message.
I had a look at the following tutorial [Setting multiple Dirichlet, Neumann, and Robin conditions β FEniCSx tutorial](https://jsdokken.com/dolfinx-tutorial/chapter3/robin_neumann_dirichlet.html) where the boundary integral vanishes for the Dirichlet part since the test function is equal to 0.
I have also consider including the boundary integraIs by using
`FacetNormal(mesh)` to compute the unit outward normal vector on the boundary and the operator `grad()` to extract the derivative in the direction normal to the boundary.
Here is the code.
```auto
n = FacetNormal(mesh)
# Weak form: Volume integrals
a_volume = (
-F_S22 * f.dx(0) * v.dx(0) * dx
- F_S11 * f.dx(1) * v.dx(1) * dx
+ 2 * F_S12 * f.dx(0) * v.dx(1) * dx
)
# Weak form: Boundary integrals
a_boundary_F22 = F_S22 * dot(grad(f), n) * v * ds # F_{,22} term
a_boundary_F11 = F_S11 * dot(grad(f), n) * v * ds # F_{,11} term
a_boundary_F12 = F_S12 * (grad(f)[0] * n[1] + grad(f)[1] * n[0]) * v * ds # F_{,12} term
# Combine volume and boundary terms
a = a_volume + a_boundary_F22 + a_boundary_F11 + a_boundary_F12
# Linear form (RHS)
L = load * v * dx
# Boundary condition
bc = DirichletBC(V, Constant(0.0), "on_boundary")
# Solve
f_solution = Function(V)
solve(a == L, f_solution, bc)
```
The results look stable.
Which is the correctly way to handle the weak form? Should I neglect the boundary integrals given the Dirichlet conditions, or should they be included? If they should be included, should anything be adjusted or included differently?
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [December 30, 2024, 9:57am UTC](https://fenicsproject.discourse.group/t/how-to-implement-first-and-second-derivatives-w-r-t-specific-directions-in-fenics/16536/9 "2024-12-30T09:57:53Z")
If you have Dirichlet conditions, any boundary integral over that boundary is ignored. See for instance
[http://jsdokken.com/FEniCS-workshop/src/form\_compilation/alternate\_form.html#lifting](http://jsdokken.com/FEniCS-workshop/src/form_compilation/alternate_form.html#lifting)