what exactly does the ufl.grad function do?

1

Hi,

what exactly does the ufl.grad function do?

I have an array A with 1076 entries.
Now I use the function:

B=ufl.grad(A)

B now has 2074 entries. Why does B have more entries than A?

2

In general, ufl.grad is a symbolic representation of the gradient of some function that depends on spatial coordinates, which is used by the form compiler to generate lower-level code for number crunching. The meaning of applying it to an array is unclear. (If I try to reproduce your MWE by interpreting it literally and passing a NumPy array to ufl.grad, I just get an “invalid type conversion” error.)

1 Like
3

Sorry, I expressed myself wrongly. I did not mean a NumPy array.

This is an MWE. That’s not part of the programme because the meshing is too big for an MWE.

import numpy as np
from dolfinx.fem import (Constant, Function, FunctionSpace, 
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, dx, ds, grad
from dolfinx import fem
from petsc4py.PETSc import ScalarType

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(mesh, ("CG", 1))
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
x = SpatialCoordinate(mesh)
g = - 4 * x[1]
f = Constant(mesh, ScalarType(-6))
L = f * v * dx - g * v * ds

dofs_R = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 1))
u_R = Function(V)
u_R.interpolate(lambda x: 2 + 2*x[1]**2)
bc_R = dirichletbc(u_R, dofs_R)
bcs = [bc_R]

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
V2 = FunctionSpace(mesh, ("CG", 2))

B=-grad(uh)
print(np.size(uh.x.array))
B2 = fem.Expression(B, V2.element.interpolation_points)
B3 = fem.Function(V2)
print(np.size(B3.x.array))

The output ist:

uh: 121 B3: 441

Shouldn’t both outputs be the same size?

4

The difference in size for the degree-of-freedom (DoF) vectors between ("CG",2) and ("CG",1) is not surprising, since a quadratic CG space has more DoFs than a linear one.

I think maybe there is some confusion here, though, since it looks like you’re trying to interpolate a vector-valued Expression into a scalar-valued space. Keep in mind that the “2” in the construction of V2 refers to the polynomial degree of the shape functions, not the dimension of functions’ values at points. This is why calling B3.interpolate(B2) results in an error. What you need is to define V2 as a VectorFunctionSpace (which can be imported from dolfinx.fem alongside FunctionSpace) to interpolate the gradient of a scalar field, which is a vector field.

Another point to keep in mind is that spatial differentiation will decrease polynomial degree and continuity. On simplices (i.e., triangles or tetrahedra), the gradient of a ("CG",1) function will be constant on each element and discontinuous between elements, i.e, a ("DG",0) (vector-valued) function. Thus, the “correct” space to use (assuming you want the exact gradient of a ("CG",1) function) is

V2 = VectorFunctionSpace(mesh, ("DG", 0))
2 Likes
5

Sorry in my code it works B3.interpolate(B2). My problem is: I need to multiply the two arrays uh and B3 together. Is there any way to align the two arrays?

6

You would then need to interpolate them into a common function space, or use a projection.