Problem with Projection

Rene_Thierry_Djoumes · June 20, 2024, 10:51am

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
    grad_u = grad(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
    grad_uh = fem.Function(V)
    problem = fem.petsc.LinearProblem(a, L)
    grad_uh = problem.solve()
    
    # Normalize the projected gradient
    norm_grad_uh = fem.Function(V)
    norm_uh_expr = grad_uh / sqrt(dot(grad_uh, grad_uh))
    
    # Interpolate the normalized gradient
    norm_grad_uh.interpolate(lambda x: norm_uh_expr.eval(x))
    
    return norm_grad_uh

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?

Thanks in advance.

francesco-ballarin · June 20, 2024, 11:25am

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.

Rene_Thierry_Djoumes · June 20, 2024, 12:45pm

Thank you for your prompt reply.
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
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
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
    grad_u = grad(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
    grad_uh = fem.Function(V)
    problem = fem.petsc.LinearProblem(a, L)
    grad_uh = problem.solve()
    
    # Normalize the projected gradient
    norm_grad_uh = fem.Function(V)
    norm_uh_expr = grad_uh / sqrt(dot(grad_uh, grad_uh))
    
    # Interpolate the normalized gradient
    #norm_grad_uh.interpolate(lambda x: norm_uh_expr.eval(x))
    norm_grad_uh.interpolate(norm_uh_expr)
    
    return norm_grad_uh


Nr = normalizeAndProj_vec(uh, domain, FiberSpace)

francesco-ballarin · June 20, 2024, 1:28pm

Write the error within “```”.

Rene_Thierry_Djoumes · June 20, 2024, 1:50pm

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
    norm_grad_uh.interpolate(norm_uh_expr)
  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()

francesco-ballarin · June 20, 2024, 3:55pm

Projection may not be needed at all.

I’ve simplified your code to

import mpi4py
import dolfinx.fem
import ufl


def normalize_gradient(u, V):
    """
    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.
    """

    grad_uh = ufl.grad(uh)
    norm_grad_u_ufl = grad_uh / ufl.sqrt(ufl.dot(grad_uh, grad_uh))

    norm_grad_u_expr = dolfinx.fem.Expression(norm_grad_u_ufl, V.element.interpolation_points())
    norm_grad_u_interp = dolfinx.fem.Function(V)
    norm_grad_u_interp.interpolate(norm_grad_u_expr)

    return norm_grad_u_interp


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, )))
norm_grad_uh_interp = normalize_gradient(uh, FiberSpace)

I haven’t tested this thoroughly (not even exporting the resulting field for visualization in paraview), but hopefully it should give you an hint on how to proceed in your case.

Rene_Thierry_Djoumes · June 21, 2024, 1:03pm

Thank very much @francesco-ballarin its works.
But I would not close the discussion now.

Topic		Replies	Views
$L^2$ Projection mathematics	5	146	February 6, 2024
Where to find 'project'-function in dolfinx? dolfinx	4	2570	November 30, 2022
Projection operator on test functions variational formulation	0	958	April 9, 2020
Find gradient of solution in dolfinx dolfinx	2	1682	June 21, 2021
Projection of values General mesh	3	108	April 15, 2024

Problem with Projection

Related topics