How to Implement First and Second Derivatives w.r.t. Specific Directions in FEniCS

[Great_Mammuth]

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:

    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

[Stein]

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)

[Great_Mammuth]

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:

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

[dokken]

@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?

Biharmonic equation: Biharmonic equation — DOLFINx 0.10.0.0 documentation 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. 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

[Great_Mammuth]

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:

   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

[dokken]

They look fine to me, but again, one avoids second derivatives in standard Galerkin problems, thus you should integrate these terms by parts.

[Great_Mammuth]

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

[Great_Mammuth]

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 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.

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?

[dokken]

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