Point source for Nedelec elements

dolfinx
1

Hi,

Can we use scifem.PointSource with Nedelec elements? At least a direct copy of the example does not work since I cannot do V.sub(1) when V is an N1curl object.

Thanks!

2

A point-source on a sub-component of a Nedelec element is not something that I’ve considered, as a Nedelec space does not have a sub-space.

The point-source would then be some vector field with a magnitude in x and y-direction.
Would this be what you would like?

Additionally, since Nedelec functions are scaled by a covariant Piola transform (Nédélec (first kind) | DefElement), How to define a finite element | DefElement.
The applied source becomes mesh dependent, as it depends on J^{-T}.

I’ve made a branch of scifem that handles this kind of point source, see:

for details

1 Like
3

I am afraid I cannot see far enough on technical implementation challenges. To the extent I get it, the point sources needed in electromagnetics should be possibly defined in a desired direction in all 2D/3D space (depending upon geometry).

As for mesh dependence of J^{-T}, I have always had to do a convergence study to trust results. As long as things stabilize for a reasonable mesh I do not see it as a problem.

Thanks a lot for your help!

4

You can test that branch of scifem, which allows you to specify a «component» for a vector-valued space such as Nedelec.

1 Like
5

I am sorry for late reply - got distracted and could not work it out right away.

I installed this new branch in my 0.9.0 installation but ran into the following error:

ValueError                                Traceback (most recent call last)
Cell In[12], line 65
     63 b = dolfinx.fem.Function(Omega)
     64 b.x.array[:] = 0
---> 65 point_source.apply_to_vector(b)

File ~/miniconda3/envs/fenics-complex-0.9/lib/python3.13/site-packages/scifem/point_source.py:210, in PointSource.apply_to_vector(self, b, recompute)
    208 for dofs, values in zip(unrolled_dofs, self._basis_values, strict=True):
    209     if isinstance(b, dolfinx.fem.Function):
--> 210         b.x.array[dofs] += values * self._magnitude
    211     elif isinstance(b, dolfinx.la.Vector):
    212         b.array[dofs] += values * self._magnitude

ValueError: operands could not be broadcast together with shapes (20,) (20,3) (20,)

This simple MWE reproduces the error

import dolfinx, ufl, gmsh, petsc4py, mpi4py, numpy as np
from dolfinx.io import gmshio
import scifem

fem_order = 2

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gmsh.clear()
gdim = 3
model_rank = 0

occ = gmsh.model.occ

fov = occ.addSphere(0, 0, 0, 1, angle3=np.pi/2)
pnt = occ.addPoint(0, 0, 0)
occ.synchronize()

frag, _ = occ.fragment([(gdim, fov)], [(0, pnt)])

# number all domains
all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")
model_rank = 0
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)

Omega = dolfinx.fem.functionspace(mesh, ("N1curl", fem_order))

# point source

geom_dtype = mesh.geometry.x.dtype
if mesh.comm.rank == 0:
    points = np.array([0, 0.0, 0], dtype=geom_dtype)
else:
    points = np.zeros((0, 3), dtype=geom_dtype)

# Next, we create the point source object and apply it to the right hand side vector.

gamma = 1.
point_source = scifem.PointSource(Omega, points, magnitude=gamma, component=0)
b = dolfinx.fem.Function(Omega)
b.x.array[:] = 0
point_source.apply_to_vector(b)
6

I just pushed a commit that makes it possible to do point-sources on 0.9 with nedelec elements.

1 Like