# Problem with Projection

Dear all,

I’m trying to project the norm of the gradient of a solution into a function space. In the FEniCs version, this operation was possible using project(). In FEniCSx, this procedure is not as simple.
I followed a discussion (Find gradient of solution in dolfinx) on this subject and created a function for projection but I always get an error.
The equation I solve is the poisson equation with the source term zero
Here’s the function I’m using and the error

def normalizeAndProj_vec(u, mesh, V):
"""
Normalizes the gradient of the solution u and projects it into the specified V.

Parameters:
u (Function): The solution whose gradient is to be normalized.
mesh (Mesh): The mesh on which the function is defined.
V (FunctionSpace): The function space into which the normalized gradient is to be projected.

Returns:
Function: The normalized gradient projected into FiberSpace.
"""

# Calculate the gradient of u

# Define the test and trial functions
v = TestFunction(V)
w = TrialFunction(V)

# Define the variational problem for the projection
a = inner(w, v) * dx
L = inner(grad_u, v) * dx

# Solve the projection problem
problem = fem.petsc.LinearProblem(a, L)



the error

  File "/home/djoumessi/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/operators.py", line 167, in inner
return Inner(a, b)
^^^^^^^^^^^
File "/home/djoumessi/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/tensoralgebra.py", line 162, in __new__
raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
ValueError: Shapes do not match: <Grad id=139755148445568> and <Argument id=139756261500464>


I was wondering if there wasn’t a simple approach to this?

The code you provided isn’t complete. The error is probably in how you define V when calling your normalizeAndProj_vec function, but we can only guess since you didn’t post the full code. See Read before posting: How do I get my question answered?

Say that V is a FE space for a scalar field, and hence u is a scalar field. When you ask grad(u), you are getting a vector field, and so you must be passing a FE space which is not V (scalar field), but the equivalent FE space of V for a vector field.

1 Like

sorry for the inconvenience.
I change the V to be a VectorFunctionSpace but I’ve another error,
Here is the MWE reproduicing the error.

Thanks

from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem import functionspace, FunctionSpace, VectorFunctionSpace
from ufl import grad, TestFunction, TrialFunction, inner, dx, sqrt, dot
# create domain and function space
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = FunctionSpace(domain, ("Lagrange", 1))
FiberSpace = VectorFunctionSpace (domain, ("Lagrange", 1))
from dolfinx import fem

uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
import numpy
# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)
import ufl
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
from dolfinx import default_scalar_type
f = fem.Constant(domain, default_scalar_type(-6))
# Bi and Linear form
L = f * v * ufl.dx

from dolfinx.fem.petsc import LinearProblem
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

# function for the projection

def normalizeAndProj_vec(u, mesh, V):
"""
Normalizes the gradient of the solution u and projects it into the specified V.

Parameters:
u (Function): The solution whose gradient is to be normalized.
mesh (Mesh): The mesh on which the function is defined.
V (FunctionSpace): The function space into which the normalized gradient is to be projected.

Returns:
Function: The normalized gradient projected into FiberSpace.
"""

# Calculate the gradient of u

# Define the test and trial functions
v = TestFunction(V)
w = TrialFunction(V)

# Define the variational problem for the projection
a = inner(w, v) * dx
L = inner(grad_u, v) * dx

# Solve the projection problem
problem = fem.petsc.LinearProblem(a, L)

Nr = normalizeAndProj_vec(uh, domain, FiberSpace)


Write the error within “”.

here is the error

WARNING:py.warnings:/home/djoumessi/Bureau/IMT4/Fenicsx_tutorial/read_mesh/minimal/MWE.py:9: DeprecationWarning: This method is deprecated. Use FunctionSpace with an element shape argument instead
FiberSpace = VectorFunctionSpace (domain, ("Lagrange", 1))

Traceback (most recent call last):
File "/home/djoumessi/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/function.py", line 397, in interpolate
_interpolate(u, cells)
File "/home/djoumessi/miniconda3/envs/fenicsx-env/lib/python3.12/functools.py", line 909, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/djoumessi/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/function.py", line 373, in _interpolate
self._cpp_object.interpolate(u, cells, nmm_interpolation_data)
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
1. (self: dolfinx.cpp.fem.Function_float64, f: numpy.ndarray[numpy.float64], cells: numpy.ndarray[numpy.int32]) -> None
2. (self: dolfinx.cpp.fem.Function_float64, u: dolfinx.cpp.fem.Function_float64, cells: numpy.ndarray[numpy.int32], nmm_interpolation_data: Tuple[List[int], List[int], List[float], List[int]]) -> None
3. (self: dolfinx.cpp.fem.Function_float64, expr: dolfinx::fem::Expression<double, double>, cells: numpy.ndarray[numpy.int32]) -> None

Invoked with: <dolfinx.cpp.fem.Function_float64 object at 0x7f5181e4ccf0>, ComponentTensor(Division(Indexed(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, quadrilateral, 1, gll_warped, unset, False), (2,)), 0), blocked element (Basix element (P, quadrilateral, 1, gll_warped, unset, False), (2,))), 3), MultiIndex((Index(8),))), Sqrt(Dot(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, quadrilateral, 1, gll_warped, unset, False), (2,)), 0), blocked element (Basix element (P, quadrilateral, 1, gll_warped, unset, False), (2,))), 3), Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, quadrilateral, 1, gll_warped, unset, False), (2,)), 0), blocked element (Basix element (P, quadrilateral, 1, gll_warped, unset, False), (2,))), 3)))), MultiIndex((Index(8),))), array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63], dtype=int32), ((), (), (), ())

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/home/djoumessi/Bureau/IMT4/Fenicsx_tutorial/read_mesh/minimal/MWE.py", line 78, in <module>
Nr = normalizeAndProj_vec(uh, domain, FiberSpace)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/djoumessi/Bureau/IMT4/Fenicsx_tutorial/read_mesh/minimal/MWE.py", line 73, in normalizeAndProj_vec
File "/home/djoumessi/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/function.py", line 402, in interpolate
self._cpp_object.interpolate(np.asarray(u(x), dtype=self.dtype), cells)
^^^^
File "/home/djoumessi/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/exproperators.py", line 340, in _call
if arg in ("+", "-"):
^^^^^^^^^^^^^^^^^
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()


Projection may not be needed at all.

import mpi4py
import dolfinx.fem
import ufl

"""
Normalizes the gradient of the solution u and interpolates it into the specified V.

Parameters:
u (Function): The solution whose gradient is to be normalized.
V (FunctionSpace): The function space into which the normalized gradient is to be interpolated.

Returns:
Function: The normalized gradient projected into FiberSpace.
"""

domain = dolfinx.mesh.create_unit_square(mpi4py.MPI.COMM_WORLD, 8, 8, dolfinx.mesh.CellType.quadrilateral)

V = dolfinx.fem.functionspace(domain, ("P", 2))
uh = dolfinx.fem.Function(V)
uh.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)

# Gradient of a continuous function may be discontinuous -> use DP instead of P
# Gradient of a polynomial of degree p is a polynomial of degree p - 1 -> decrease the polynomial order
FiberSpace = dolfinx.fem.functionspace(domain, ("DP", 1, (domain.topology.dim, )))
`