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.dot(ufl.grad(u), ufl.grad(u))))
But it’s wrong, I got the errors
File "/mnt/g/workdir/MKM/NO_mkm/fenics/test_refine2.py", line 16, in <module>
F = dolfinx.fem.form(ufl.sqrt(ufl.dot(ufl.grad(u), ufl.grad(u))))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/mnt/d/software_install/fenics/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 188, in form
return _create_form(form)
^^^^^^^^^^^^^^^^^^
File "/mnt/d/software_install/fenics/lib/python3.12/site-packages/dolfinx/fem/forms.py", 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/expr.py", line 372, in __iter__
for i in range(len(self)):
^^^^^^^^^
File "/mnt/d/software_install/fenics/lib/python3.12/site-packages/ufl/core/expr.py", line 368, in __len__
raise NotImplementedError("Cannot take length of non-vector expression.")
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)
u.interpolate(func)
F = dolfinx.fem.form(ufl.dot(ufl.dot(ufl.grad(u), 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 ufl.dot, and you are missing an integration domain.
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.
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(
MPI.COMM_WORLD,
points=((-3.0, -1.0), (3.0, 1.0)),
n=(30, 10),
cell_type=dolfinx.mesh.CellType.triangle)
#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)
u.interpolate(func)
#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.dot(ufl.grad(u), ufl.grad(u))),
Vg.element.interpolation_points())
grad_norm = dolfinx.fem.Function(Vg)
grad_norm.interpolate(F)
#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
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
warped = grid.warp_by_scalar()
plotter.add_mesh(warped)
plotter.show()
