How to get the first order derivative of the solution and store it in a numpy array

I/O
1

Dear All, my program works fine, but my question is how to get the first order derivative of the solution and store it in a numpy array. The core parts of my program is as follows


...
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(self.bl, self.br), nx=self.num)
...
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = FunctionSpace(msh, ufl.MixedElement(P1, P1, P1))
...
# Define functions 
u = Function(V)
 # Split system functions to access components
u1, u2, u3 = split(u)
...
V0, dofs0 = V.sub(0).collapse()
V1, dofs1 = V.sub(1).collapse()
V2, dofs2 = V.sub(2).collapse()
...
# # Create nonlinear problem and Newton solver 
problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
...
out = solver.solve(u)
sols[:,1] = u.x.array[dofs0]
sols[:,2] = u.x.array[dofs1]
sols[:,3] = u.x.array[dofs2]

Up to this point, I have obtained the solutions in the form of arrays (sols), how do I get the first order derivative of the solution next? I had tried the following methods, like this

udx1= u.dx(0).x.array[dofs0]

but it do not work.

Any suggestions would be greatly appreciated!

2

Use dolfinx.fem.Expression to interpolate into an appropriate space, and then access the dolfinx.fem.Function.x.array.

See Solving a time-dependent problem — FEniCS 22 tutorial

Note that Getting the derivatives at the degrees of freedom is not well defined.

3

Dear dokken

Thank you very much for your thread, with my current understanding of fenicsx I am not able to solve my problem through the tutorials, the problem in the tutorials has only one function to be solved, whereas mine has three, so they are significantly different.

Can you give me some specific actionable advices based on my problem above?

4

Only if you provide a minimal reproducible example. I cannot run the code above as it is lacking tons of definitions, and would require a fair lot of work and assumptions from my end to Get to something working.

5

Thank you dokken, following is my example,



'''
coordinate:
     x->[0,10]
equations:
    -u1'' - 4* u2'*u1' + 2*a*u1- 8*u1**3 + 6*u1**5 =0 , where the ' = d/dx,
   # -u2''- 2/3 *u1'**2 =0 splited into
     u3' +  2/3 *u1'**2 =0
     u3  =  u2'
boundary conditions
    u1(0) =0, # u1'(10) =0
    u2(0) =0, #u2'(0) =0
    u3(0) =0


'''

import numpy as np

import ufl
from dolfinx import mesh, fem,log
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, split, TestFunctions

from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
# parmeters
r_0 = 0.0
r_b = 10.0  # the coordinate of right boundary


a=0.2
b=2.0
c=1.0
beta=0.2
num =1000
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(r_0, r_b), nx=num-1)

fdim = msh.topology.dim - 1
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = FunctionSpace(msh, ufl.MixedElement(P1, P1,P1))

# the left boundary
facets_l = mesh.locate_entities_boundary(msh, dim=fdim, marker=lambda x: np.isclose(x[0], r_0))

dofs_u1 = fem.locate_dofs_topological(V=V.sub(0), entity_dim=fdim, entities=facets_l)
dofs_u2 = fem.locate_dofs_topological(V=V.sub(1), entity_dim=fdim, entities=facets_l)
dofs_u3 = fem.locate_dofs_topological(V=V.sub(2), entity_dim=fdim, entities=facets_l)

# the Dirichleb boundary conditions
bc_u1 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u1, V=V.sub(0))  # u1 = phi
bc_u2 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u2, V=V.sub(1))  # u2 = A
bc_u3 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u3, V=V.sub(2))  # u2 = A

bcs = [bc_u1, bc_u2, bc_u3]

# Define functions
u = Function(V)
# Split system functions to access components
#phi, m, sigma = u1, u2, u3
u1, u2,u3 = split(u)

# Define test functions
v1, v2, v3 = TestFunctions(V)

# Define variational problem

F1 = v1.dx(0)*u1.dx(0) + (-4.0 * u2.dx(0)*u1.dx(0) +2.0*a*u1- 4.0*b*u1**3 + 6.0*c*u1**5)*v1
F2 = (u3.dx(0)+ 2.0/3.0*u1.dx(0)**2)*v2
F3 = (u3-u2.dx(0))*v3
F = (F1+F2+F3)*dx
# guess solutions
x=np.linspace(r_0,r_b,num)
y1=[0,1,1,1,1,1]
y3=[0,0,-0.8,-1.2,-1.6,-2.0]
x1=np.linspace(r_0,r_b,num=6)
f1=interp1d(x1,y1,kind='linear')
f3=interp1d(x1,y3,kind='linear')
y1_ = f1(x)
y2_ = np.diff(y1_)
y2_ = np.insert(y2_,0,0)
y3_ = f3(x)
#yguess=guess(x)
V0, dofs0 = V.sub(0).collapse()
V1, dofs1 = V.sub(1).collapse()
V2, dofs2 = V.sub(2).collapse()
u.x.array[dofs0]=y1_ # yguess[0]
u.x.array[dofs1]=y2_ # yguess[1]
u.x.array[dofs2]=y3_ # yguess[2]
# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-8
solver.nonzero_initial_guess = True
solver.max_it = 100
solver.report = True

out = solver.solve(u)
u0_ = u.x.array[dofs0]
u1_ = u.x.array[dofs1]
u2_ = u.x.array[dofs2]

plt.figure(2,dpi=300)
plt.plot (x, u0_,label="u1")
plt.plot (x, u1_,label="u2")
plt.legend()
plt.show()

plt.close ( )

My questions are
(1) how to get the first order derivative of the u0_?
(2) how by using your advice " Use dolfinx.fem.Expression to interpolate into an appropriate space, and then access the dolfinx.fem.Function.x.array." to improve our codes?

Thanks!

6

@dokken Dear Dokken, I have given an example above, can you give me some suggestions again?

7

Hello @dokken, I would like to know whether your last statement is true when we use 2nd order or higher CG elements. Why would the derivative be faulty in this case?
I need to compute accurately an expression involving a first order derivative and I suspect I wasn’t able to do so (as can be seen in my recent post history).

8

The derivatives are faulty if you evaluate them at a boundary between elements (facet or edges). They are only well defined on the boundary of a given cell, i.e say cell 0 and cell 1 share a point (0,0). Then the derivative is cell defined at (0,0) for both cell 0 and 1, but they might give different results.

1 Like
9

I see, thanks a lot.
ChatGPT suggests to compute the derivative at ‘‘quadrature points’’ within each cell, rather than at the vertices or edges. It says those are the points where integration is performed and that the derivatives are well defined at those points. Does it make any sense?
If so, is there a way to implement this with Fenicsx?

10

@raboynics @dokken Dear Raboynics and Dokken, I am confused with your above comments.

11

Yes, that makes sense.

See for instance: Numerical values from ufl.SpatialCoordinate - #2 by dokken

1 Like
12

I think I have found out the method to solve my question by the following codes

expr = fem.Expression(grad(u1), V1.element.interpolation_points())
w = Function(V)
u1_dx = w.sub(0).collapse()
u1_dx.interpolate(expr)

The whole example codes are be changed as follows:

'''
coordinate:
     x->[0,10]
equations:
    -u1'' - 4* u2'*u1' + 2*a*u1- 8*u1**3 + 6*u1**5 =0 , where the ' = d/dx,
   # -u2''- 2/3 *u1'**2 =0 splited into
     u3' +  2/3 *u1'**2 =0
     u3  =  u2'
boundary conditions
    u1(0) =0, # u1'(10) =0
    u2(0) =0, #u2'(0) =0
    u3(0) =0


'''

import numpy as np

import ufl
from dolfinx import mesh, fem,log
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, split, TestFunctions,grad

from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
# parmeters
r_0 = 0.0
r_b = 10.0  # the coordinate of right boundary


a=0.2
b=2.0
c=1.0
beta=0.2
num =1000
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(r_0, r_b), nx=num-1)

fdim = msh.topology.dim - 1
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = FunctionSpace(msh, ufl.MixedElement(P1, P1,P1))

# the left boundary
facets_l = mesh.locate_entities_boundary(msh, dim=fdim, marker=lambda x: np.isclose(x[0], r_0))

dofs_u1 = fem.locate_dofs_topological(V=V.sub(0), entity_dim=fdim, entities=facets_l)
dofs_u2 = fem.locate_dofs_topological(V=V.sub(1), entity_dim=fdim, entities=facets_l)
dofs_u3 = fem.locate_dofs_topological(V=V.sub(2), entity_dim=fdim, entities=facets_l)

# the Dirichleb boundary conditions
bc_u1 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u1, V=V.sub(0))  # u1 = phi
bc_u2 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u2, V=V.sub(1))  # u2 = A
bc_u3 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u3, V=V.sub(2))  # u2 = A

bcs = [bc_u1, bc_u2, bc_u3]

# Define functions
u = Function(V)
# Split system functions to access components
#phi, m, sigma = u1, u2, u3
u1, u2,u3 = split(u)

# Define test functions
v1, v2, v3 = TestFunctions(V)

# Define variational problem

F1 = v1.dx(0)*u1.dx(0) + (-4.0 * u2.dx(0)*u1.dx(0) +2.0*a*u1- 4.0*b*u1**3 + 6.0*c*u1**5)*v1
F2 = (u3.dx(0)+ 2.0/3.0*u1.dx(0)**2)*v2
F3 = (u3-u2.dx(0))*v3
F = (F1+F2+F3)*dx
# guess solutions
x=np.linspace(r_0,r_b,num)
y1=[0,1,1,1,1,1]
y3=[0,0,-0.8,-1.2,-1.6,-2.0]
x1=np.linspace(r_0,r_b,num=6)
f1=interp1d(x1,y1,kind='linear')
f3=interp1d(x1,y3,kind='linear')
y1_ = f1(x)
y2_ = np.diff(y1_)
y2_ = np.insert(y2_,0,0)
y3_ = f3(x)
#yguess=guess(x)
V1, dofs0 = V.sub(0).collapse()
V2, dofs1 = V.sub(1).collapse()
V3, dofs2 = V.sub(2).collapse()
u.x.array[dofs0]=y1_ # yguess[0]
u.x.array[dofs1]=y2_ # yguess[1]
u.x.array[dofs2]=y3_ # yguess[2]
# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-8
solver.nonzero_initial_guess = True
solver.max_it = 100
solver.report = True

out = solver.solve(u)
u1_ = u.x.array[dofs0]
u2_ = u.x.array[dofs1]
u3_ = u.x.array[dofs2]

# u1_ = u.sub(0).collapse()

expr = fem.Expression(grad(u1), V1.element.interpolation_points())
w = Function(V)
u1_dx = w.sub(0).collapse()
u1_dx.interpolate(expr)


expr2 = fem.Expression(grad(u2), V2.element.interpolation_points())

u2_dx = w.sub(1).collapse()
u2_dx.interpolate(expr2)

plt.figure(1)
plt.plot (x, u1_,label="u1")
plt.plot (x, u2_,label="u2")
plt.plot (x, u3_,label="u3") # the solution of u3 is the u2_dx theoritical
plt.plot (x, (u3_ - u2_dx.x.array[:]), label="u3-u2_dx")
plt.plot (x, u1_dx.x.array[:],label="u1_dx")
plt.plot (x, u2_dx.x.array[:],label="u2_dx")
plt.legend()
plt.show()

plt.close ( )

The output fig is
Figure_1

13

I would appreciate if @dokken would comment. This seems to be the usual way I had done this (see for example Nedelec elements or Lagrange? What about the order? - #2 by raboynics) however this lead me to strange results (result seems to converge with mesh refinement but not with higher polynomial degrees).

This way of computing derivatives seems easier to implement than the one you just linked to, dokken (Numerical values from ufl.SpatialCoordinate - #2 by dokken), where it is question to define a quadrature scheme with basix.

14

As I’ve already stated: Evaluating the derivative of a Lagrange finite element function at the element boundary is not well defined. For P1, I would do the following:

import matplotlib.pyplot as plt
from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
u = dolfinx.fem.Function(V)

u.interpolate(lambda x: np.sin(3*x[0]))


Q = dolfinx.fem.functionspace(mesh, ("DG", 0))
q = dolfinx.fem.Function(Q)
expr = ufl.grad(u)[0]
q.interpolate(dolfinx.fem.Expression(expr, Q.element.interpolation_points()))


x = V.tabulate_dof_coordinates()
x_Q = Q.tabulate_dof_coordinates()
x_sort = np.argsort(x[:, 0])
plt.plot(x[x_sort, 0], u.x.array[x_sort], "-ro", label="u")
plt.plot(x_Q[:, 0], q.x.array, "o", label="q")
plt.legend()
plt.savefig("figure.png")

As a DG 0 element has its degree of freedom at the midpoint (and the value is constant per element).
figure

1 Like
15

Thanks a lot dokken. I think you made it clear that you correctly evaluate u and grad(u) using your way, where u is a piece-wise linear function defined (or evaluated) at some dofs, whereas grad(u) is its derivative and is piece-wise constant, and it is evaluate at the mid points between the dofs coordinates of u.

However, it seems to me that zhaozhh correctly computed his derivatives in his last post. He interpolates this expression: expr = fem.Expression(grad(u1), V1.element.interpolation_points()) into a suitable space, and you do something very similar: q.interpolate(dolfinx.fem.Expression(expr, Q.element.interpolation_points())).
He got it right, right?

16

He used

and V1 is not a suitable space, as it is degrees of freedom are at vertices, and the interpolation is not well defined.

1 Like
17

Thanks again a lot dokken, now I clearly see what’s wrong. I didn’t realize picking a degree 0 space would have the dofs coordinates at different points than a degree 1 space.
This leaves me with a question. If I pick a degree n for u and a degree n-1 for grad(u), do I ensure that grad(u) is well defined? I guess not, but I am not 100 percent sure.

18

The gradient is always well defined within a given cell when using a DG space of degree n-1 for the gradient. The problem is that when you evaluate at a point that are in different cells, you get different answers depending on what cell you choose. You can for instance evaluate the function in all cells with a point and compute an average to get a well defined function for DG spaces with dofs at vertices, edges and facets.

1 Like
19

Hi @dokken,

Thanks a lot for the detailed explaination on how to compute the derivatives at the center of the elements.
What would you recommend for cases in where I would in fact need to derivative values at the DOF points of a Lagrange-1 element? (i.e, I need to have the gradients at the vertices of a 2D mesh)

Could it be possible to project or interpolate the DG-0 derivative results into a Lagrange-1 function?

20

Hi again,

I wonder if interpolating the derivative function (DG -0) back into a function from the primary function space (Lagrange,1) makes any sense?

I have something like this:

Q = functionspace(domain, ("DG",0))
m1_dx = fem.Function(Q)

m1_dx.interpolate(fem.Expression(ufl.grad(m1)[0], Q.element.interpolation_points()))
##We then try to interpolate back into the Lagrange 1 function space
m1_dx_l1 = fem.Function(V)

m1_dx_l1.interpolate(m1_dx)

In this case V is a Lagrange-1 function space and m1 is a function from this function space.

Would this approach make sense? Can I trust the results from interpolating back from the DG-0 to the Lagrange-1 function space?