Calculte the norm of a gradient

Hello. I want to calculate the \sqrt{\nabla u\cdot \nabla u} of a function u, I write the variational form as the following

F = dolfinx.fem.form(ufl.sqrt(, ufl.grad(u))))

But it’s wrong, I got the errors

File "/mnt/g/workdir/MKM/NO_mkm/fenics/", line 16, in <module>
    F = dolfinx.fem.form(ufl.sqrt(, ufl.grad(u))))
  File "/mnt/d/software_install/fenics/lib/python3.12/site-packages/dolfinx/fem/", line 188, in form
    return _create_form(form)
  File "/mnt/d/software_install/fenics/lib/python3.12/site-packages/dolfinx/fem/", line 185, in _create_form
    return list(map(lambda sub_form: _create_form(sub_form), form))
  File "/mnt/d/software_install/fenics/lib/python3.12/site-packages/ufl/core/", line 372, in __iter__
    for i in range(len(self)):
  File "/mnt/d/software_install/fenics/lib/python3.12/site-packages/ufl/core/", line 368, in __len__
    raise NotImplementedError("Cannot take length of non-vector expression.")

How should I write it ?

I suspect u is a vector quantity. Then grad(u) is a tensor, and the dot product is not appropriate. Use ufl.inner instead,

Sorry for the confusion, in fact u is a scaler function. Here is my complete codes. I want to locate the cells on which \parallel \nabla u\parallel is large and refine the grids on the cells.

# -*- coding:utf-8 -*-
import ufl 
import dolfinx
from mpi4py import MPI
import numpy as np 

domain = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, points=((-3.0, -1.0), (3.0, 1.0)), n=(30, 10), cell_type=dolfinx.mesh.CellType.triangle)

def func(x):
    return np.sin(2*(x[0]+x[1])) * np.exp(-x[0]**2)

V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1))
u = dolfinx.fem.Function(V)

F = dolfinx.fem.form(, ufl.grad(u))))
grad_norm = dolfinx.fem.assemble_vector(F)

I expect that I can get a vector of \parallel \nabla u\parallel “grad_norm”, and then refine the grids…

A norm returns a scalar value.
Thus you should use dolfinx.fem.assemble_scalar
Please also note that your form no longer makes sense, there are too many, and you are missing an integration domain.

I think I messed up some conceptions, such as “form” and “expression”.
I tried this code and it works

F = dolfinx.fem.Expression(ufl.sqrt(, ufl.grad(u))), V.element.interpolation_points())
grad_norm = dolfinx.fem.Function(V)

Hello, dokken. Can I do this to a function (e.g. the grad_norm I calculated above) : locating the 30% of the cells with the highest grad_norm values and get their local (or global) index, so that I can use the dolfinx.mesh.compute_incident_entities function to get the edges that I should refine ?

I would probably project the grad norm into a DG-0 space, giving you the integrated norm per cell.

You could interpolate into a DG-0 space (evaluating the norm of the gradient at the midpoint of a cell).

I am suggesting these as the gradient norm is not well defined at vertices (it has different values in different cells).

With either of the approaches above, you can then access the underlying array, Which has a 1-1 correspondence the the cell indices.
You can use numpy to find the index of the N largest values, and finally use compute_incident_entities to get the edges.

Ok, thank you a lot. I will try these suggestions. :handshake:

Here is an example code, calculates the norm of gradient and refines the 30% of cells with the highest gradient norm.

# -*- coding:utf-8 -*-
import ufl 
import dolfinx
from mpi4py import MPI
import numpy as np 
import pyvista
from dolfinx import plot

domain = dolfinx.mesh.create_rectangle(
    points=((-3.0, -1.0), (3.0, 1.0)), 
    n=(30, 10), 

#we use this function as an example
def func(x):
    return np.sin(2*(x[0]+x[1])) * np.exp(-x[0]**2)

V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1))
u = dolfinx.fem.Function(V)

#calculate the norm of gradient. Note that, it is a function instead of form
#When the 1st order Lagrange finite element is used, u is a piecewise linear function, gradient 
#will be discontinuous at the vertex. Here, we use zero order Discontinuous Galerkin (DG) finite element
#There is a DOF at the center of the cell, it has a one-to-one correspondence with the cell.
Vg = dolfinx.fem.FunctionSpace(domain, ("DG", 0))
F = dolfinx.fem.Expression(
    ufl.sqrt(, ufl.grad(u))), 
grad_norm = dolfinx.fem.Function(Vg)

#sort the grad_norm and get the cell index
order = np.argsort(grad_norm.x.array)
cell_index = order[-int(0.3*order.size):-1]

#cell_index now contains the cell which we should refine
#now we get the edges of these cells
domain.topology.create_connectivity(1, 2)
edge_index = dolfinx.mesh.compute_incident_entities(domain.topology, cell_index.astype(np.int32), 2, 1)

#call refine
new_domain = dolfinx.mesh.refine(domain, edge_index, True)
new_V = dolfinx.fem.FunctionSpace(new_domain, ("Lagrange", 1))
new_u = dolfinx.fem.Function(new_V)

#Finally, we interpolate the old function onto the new mesh
interp_data= dolfinx.fem.create_nonmatching_meshes_interpolation_data(
    new_V.mesh._cpp_object, new_V.element, V.mesh._cpp_object, 1.0e-8)

new_u.interpolate(u, nmm_interpolation_data=interp_data)

#you can plot the function. As you can see, the cells on which gradient norms are large are refined. 
cells, types, x = plot.vtk_mesh(new_V)
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["u"] = new_u.x.array.real
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
warped = grid.warp_by_scalar()