Derivative in dolfinx

Hi, I solved an electrostatic problem and exported the electric potential. The electric potential is smooth and then its numeric derivative was calculated to get the electric field. The problem is that the electric field is not smooth now. So I think if it would be better to directly calculate the derivative in dolfinx. How to calculate derivative of a function u on domain V with regard to z in dolfinx for example?

Not everything can make a smallest example. I give up. But how to calculate derivative in dolfinx?

I don’t understand what derivative you want to calculate.
Is it a spatial derivative of a `dolfinx.fem.Function`?
Is it a derivative of an ufl-Expression with respect to a variable?

If `u` is a function in say a first order lagrange space, then its derivative in `z`-direction is a piecewise continuous function on every cell, thus it is no longer continuous.
To get such derivatives, you can use `u.dx(i)` and interpolate with `dolfinx.fem.Expression` into a suitable function space (“DG” 0 in the case of `u` being in P1).

Dear dokken,
Thanks for your help. The solution u is a second order lagrange function defined as

``````V = FunctionSpace(domain, ("CG", 3))
u = TrialFunction(V)
``````

Then the electric field Ez along z axis I calculated as:

``````W = FunctionSpace(domain, ("DG", 2))

expr = Expression(uh.dx(2), W.element.interpolation_points())

w = Function(W)

w.interpolate(expr)

tol = 0.001 # Avoid hitting the outside of the domain
z = np.linspace(0 + tol, 31.6 - tol, 101)
points = np.zeros((3, 101))
points[2] = z
u_values = []
bb_tree = geometry.BoundingBoxTree(domain, domain.topology.dim)

cells = []
points_on_proc = []
cell_candidates = geometry.compute_collisions(bb_tree, points.T)
colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points.T)
for i, point in enumerate(points.T):
points_on_proc.append(point)

points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_values = w.eval(points_on_proc, cells)
``````

Is it correct?

Yes, this would be one way of doing it. Note that if a point is on a boundary between two cells, there will be some ambiguity as to which value it should take (That is decided by whatever cell is first in `colliding_cells` for the given point).

Thanks again. The derivative, that is the electric field, is much smooth now.

Dear dokken, the electric field calculated from the derivative of the potential is still not quite smooth.

Please make a minimal reproducible example that illustrates what you get and what you expect.

The file used to calculate electric field:

``````from dolfinx.io import gmshio, XDMFFile

from dolfinx.fem import (Function, FunctionSpace, locate_dofs_topological,

dirichletbc, Constant, Expression, form, assemble_scalar, assemble_vector)

from dolfinx.fem.petsc import LinearProblem

from dolfinx import geometry

from dolfinx.mesh import compute_incident_entities

from mpi4py import MPI

from ufl import TrialFunction, TestFunction, dot, grad, dx, Measure

from petsc4py.PETSc import ScalarType

from scipy.special import i0, iv

import numpy as np

x0 = 3.4186

y0 = 0.75

alpha = np.arctan(y0 / x0)

beta = np.arctan(x0 / y0)

voltage = 0.13395 / 2

l = 3.

m = 2.4

def out_boundary(x):

theta = np.arctan(x[1] / x[0])

return -voltage + (theta - alpha) / (beta - alpha) * voltage * 2

mesh_comm = MPI.COMM_WORLD

gdim = 3

gmsh_model_rank = 0

domain, cell_markers, facet_markers = gmshio.read_from_msh('cell.msh', mesh_comm, rank=gmsh_model_rank, gdim=gdim)

V = FunctionSpace(domain, ("CG", 3))

u = TrialFunction(V)

v = TestFunction(V)

L = Constant(domain, ScalarType(0)) * v * dx

print(np.unique(facet_markers.values))

pos_pole_facets = facet_markers.indices[facet_markers.values == 32]

pos_pole_bdofs = locate_dofs_topological(V, 2, pos_pole_facets)

neg_pole_facets = facet_markers.indices[facet_markers.values == 33]

neg_pole_bdofs = locate_dofs_topological(V, 2, neg_pole_facets)

out_facets = facet_markers.indices[facet_markers.values == 34]

out_bdofs = locate_dofs_topological(V, 2, out_facets)

int_facets = facet_markers.indices[facet_markers.values == 35]

int_bdofs = locate_dofs_topological(V, 2, int_facets)

u_O = Function(V)

boundary_cells = compute_incident_entities(domain, out_facets, domain.topology.dim-1, domain.topology.dim)

u_O.interpolate(out_boundary, boundary_cells)

bc_O = dirichletbc(u_O, out_bdofs)

bc_N = dirichletbc(ScalarType(-voltage), neg_pole_bdofs, V)

bc_P = dirichletbc(ScalarType(voltage), pos_pole_bdofs, V)

bcs = [bc_O, bc_N, bc_P]

problem = LinearProblem(a, L, bcs=bcs)

uh = problem.solve()

W = FunctionSpace(domain, ("DG", 2))

expr = Expression(uh.dx(2), W.element.interpolation_points())

w = Function(W)

w.interpolate(expr)

tol = 0.001 # Avoid hitting the outside of the domain

z = np.linspace(0 + tol, l - tol, 1001)

points = np.zeros((3, 1001))

points[2] = z

u_values = []

bb_tree = geometry.BoundingBoxTree(domain, domain.topology.dim)

cells = []

points_on_proc = []

cell_candidates = geometry.compute_collisions(bb_tree, points.T)

colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points.T)

for i, point in enumerate(points.T):

points_on_proc.append(point)