# Evaluation of grad(u) at a point x\_ref for variational formulation

**URL:** <https://fenicsproject.discourse.group/t/evaluation-of-grad-u-at-a-point-x-ref-for-variational-formulation/16020>\
**Category:** variational formulation\
**Tags:** linear-algebra\
**Created:** [October 10, 2024, 7:40pm UTC](https://fenicsproject.discourse.group/t/evaluation-of-grad-u-at-a-point-x-ref-for-variational-formulation/16020 "2024-10-10T19:40:17Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Hell](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/hell/32/7664_2.png) [@Hell](https://fenicsproject.discourse.group/u/Hell)\
**Post date:** [October 10, 2024, 7:40pm UTC](https://fenicsproject.discourse.group/t/evaluation-of-grad-u-at-a-point-x-ref-for-variational-formulation/16020/1 "2024-10-10T19:40:17Z")

</div>

Hello everyone,

I am trying to solve the following problem in its variational form:

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/f33ee072c46e7c548565edad1ceeb30d69e7fc6c.png)  
In the first step, I want to solve the homogeneous eigenvalue problem, where the right-hand side is set to 0. I am using SLEPc for this. In the second step, I plan to take the computed eigenvalue, \omega\_0, and substitute it into the right-hand side to solve the inhomogeneous equation.

Now, to my issue: I am unsure how to evaluate \nabla \hat{p} at (x\_{\text{ref}}) in the variational form. Is there a way to implement something like a delta function that equals 1 at x\_{\text{ref}} and 0 elsewhere?

I also have another problem related to the shapes in the matrix definition for matrix D:

```auto
File "/home/hell/anaconda3/envs/env3-10-complex/lib/python3.10/site-packages/ufl/tensoralgebra.py", line 166, in __new__
    raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
ValueError: Shapes do not match: <ComponentTensor id=139710411729280> and <Argument id=139710413172800>

```

I’m working with FEniCSx 0.8  
Thank you in advance. If something is unclear, I could provide more information.

Code:

```auto
from dolfinx.fem import locate_dofs_geometrical, functionspace, Function, dirichletbc
from mpi4py import MPI
from numpy.array_api import complex128
from petsc4py import PETSc
from ufl import TrialFunction, TestFunction, dx, inner, grad, SpatialCoordinate
from dolfinx import default_scalar_type, fem
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from dolfinx.mesh import CellType, create_rectangle, locate_entities
import numpy as np

# ---------- Domain definition ----------
w = 0.5
h = 0.1
nx = 100
ny = 10

cell_w = w / nx

msh = create_rectangle(
    MPI.COMM_WORLD, np.array([[0, 0], [w, h]]), np.array([nx, ny]), CellType.quadrilateral
)

cell_w = w / nx

msh = create_rectangle(
    MPI.COMM_WORLD, np.array([[0, 0], [w, h]]), np.array([nx, ny]), CellType.quadrilateral
)

# ---------- Define functions and area for c1, c2, rho1, rho2 and flame region ----------
Q = functionspace(msh, ("DG", 0))

def Omega_0(x):
    return x[0] <= w/2

def Omega_1(x):
    return x[0] >= w/2

def flame_region(x):
    return np.logical_and(x[0] >= w/2 - cell_w, x[0] <= w/2 + cell_w)

def ref_point(x):
    return np.logical_and(x[0] >= w/2-cell_w, x[0] <= w/2)

cells_1 = locate_entities(msh, msh.topology.dim, Omega_0)
cells_2 = locate_entities(msh, msh.topology.dim, Omega_1)
cells_flame_region = locate_entities(msh, msh.topology.dim, flame_region)
ref_point_region = locate_entities(msh,msh.topology.dim, ref_point)

c = Function(Q)
c.x.array[cells_1] = np.full_like(cells_1, 347.18, dtype=default_scalar_type)
c.x.array[cells_2] = np.full_like(cells_2, 694.36, dtype=default_scalar_type)

rho = Function(Q)
rho.x.array[cells_1] = np.full_like(cells_1, 4, dtype=default_scalar_type)
rho.x.array[cells_2] = np.full_like(cells_2, 1, dtype=default_scalar_type)

# ---------- Define functionspace, functions and BCs ----------
degree = 1
V = functionspace(msh, ('Lagrange', degree))
u = TrialFunction(V)
v = TestFunction(V)
e = Function(V)
o = Function(V)

x = SpatialCoordinate(msh)
dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], w))
bc = [dirichletbc(default_scalar_type(0), dofs, V)]

# Eigenvalues from SLEPc
omega0 = [(1709.4992392760605+3.1592511086743223e-17j), (4363.241033918917-3.6488298900508794e-17j), (7018.447703170593-2.4364863230542956e-15j), (10443.252332619693-3.7277152271679734e-15j)]
rho_ref = 4
n_ = 0.01
delta_f = cell_w
tau = 10e-4

k = 1 / rho_ref * n_ / delta_f * np.exp(1j * tau * omega0[0])
e.x.array[cells_flame_region] = np.full_like(cells_flame_region, k, dtype=complex128)
o.x.array[ref_point_region] = np.full_like(ref_point_region,k, dtype=complex128)
# ---------- Define D matrix for flame model ----------
A = inner((1 / rho) * grad(u), grad(v)) * dx
M = inner(1 / (c ** 2 * rho) * u, v) * dx
A = assemble_matrix(fem.form(A), bcs=bc)
A.assemble()
M = assemble_matrix(fem.form(M), bcs=bc)
M.assemble()
D = inner(grad(u)*e, u) * dx
D = assemble_vector(fem.form(D))
D.assemble()
P = A + D

# ---------- Solve eigenvalue problem with SLEPc ----------

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 10, 2024, 8:40pm UTC](https://fenicsproject.discourse.group/t/evaluation-of-grad-u-at-a-point-x-ref-for-variational-formulation/16020/2 "2024-10-10T20:40:38Z")

</div>

There are various topics on this, see:

> [@How to construct a sparse matrix used for calculating function values at multiple points](https://fenicsproject.discourse.group/t/how-to-construct-a-sparse-matrix-used-for-calculating-function-values-at-multiple-points/13081/4):
>
> What you are searching for is how to create a point source. Following is code for getting the basis functions evaluated at a point in the cell: # Author: Jørgen S. Dokken # # Create a point source import dolfinx from mpi4py import MPI import numpy as np import ufl def compute\_cell\_contributions(V, points): # Determine what process owns a point and what cells it lies within \_, \_, owning\_points, cells = dolfinx.cpp.geometry.determine\_point\_ownership( mesh.\_cpp\_object, points, …

and

> [@Eval in Expression freezes in parallel](https://fenicsproject.discourse.group/t/eval-in-expression-freezes-in-parallel/13840/2):
>
> It freezes as you are not sending in the appropriate MPI communicator to the expression. Since you are not sending in a communicator (which in your case should be MPI.COMM\_SELF, it is using the mesh communicator, thus waiting for a sync step between all processes that never happens. See: [dolfinx.fem — DOLFINx 0.7.3 documentation](https://docs.fenicsproject.org/dolfinx/v0.7.3/python/generated/dolfinx.fem.html?highlight=expression#dolfinx.fem.Expression) i.e. E = Expression(expr, x\_ref, comm=MPI.COMM\_SELF) does the trick

I think this is also implemented in helmholtx by @Ekrem_Ekici : [helmholtz-x/helmholtz\_x/flame\_matrices.py at 3a20260bb51a48a8c315f6c02277465ddd729e01 · ekremekc/helmholtz-x · GitHub](https://github.com/ekremekc/helmholtz-x/blob/3a20260bb51a48a8c315f6c02277465ddd729e01/helmholtz_x/flame_matrices.py#L129-L160)  
which I guess is what you want.

---

<div class="post-metadata">

**Author:** ![Hell](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/hell/32/7664_2.png) [@Hell](https://fenicsproject.discourse.group/u/Hell)\
**Post date:** [October 11, 2024, 6:46am UTC](https://fenicsproject.discourse.group/t/evaluation-of-grad-u-at-a-point-x-ref-for-variational-formulation/16020/3 "2024-10-11T06:46:32Z")

</div>

Thank you very much.  
Helmholtz-x looks very promising for my use.