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?
Please supply a minimal reproducible example that illustrates your issue
Not everything can make a smallest example. I give up. But how to calculate derivative in dolfinx?
Why cannot upload file?
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).
See for instance:
https://jsdokken.com/fenics22-tutorial/heat_eq.html#post-processing-without-projections
and
https://fenics.readthedocs.io/projects/ufl/en/latest/manual/form_language.html?highlight=dx#basic-spatial-derivatives
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):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
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.
My mesh file can be downloaded from https://musetransfer.com/s/h28z6rsdi
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)
a = dot(grad(u), grad(v)) * dx
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):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_values = w.eval(points_on_proc, cells)
np.savetxt('Ez.txt', np.c_[ points[2], u_values])
This is an electrostatic problem, and the potential is first solved. Then, calculate the electric field along the z-axis by making the derivative of the electric potential relative to z. I hope the electric field saved in the file Ez.txt is smooth, but when zooming out, it is not too smooth. The first column of Ez.txt in the figure is x, and the second column is y.
Hello. I don’t understand why you use the “scalar discontinuous Lagrange finite elements” (DG) to interpolate the derivative. Will the discontinuity disappear if you change it to “CG” ?