# Unable to compute the curl of a vector field

**URL:** <https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979>\
**Category:** dolfinx\
**Created:** [June 13, 2024, 8:38pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979 "2024-06-13T20:38:30Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 13, 2024, 8:38pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/1 "2024-06-13T20:38:30Z")

</div>

Hi, I work on a 2D problem and my use-case code is very long, I know I must create a MWE but I greatly struggle, as I am not used to employ a mesh that doesn’t come from Gmsh. Currently my attempt is unfinished and I will post it at the end (so maybe we can finish the MWE and I might use it for my next questions).

My problem lies in errors when computing the curl of a vector field. The relevant part of my code causing me problem is:

```auto
(...)
deg = 4
# Define ME function space
el = FiniteElement("CG", mesh.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(mesh, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)
(...)
S_xx = 3e-3
S_yy = S_xx * 10
rho = 1
sigma = 1.0 / rho
κ = 1.8
Seebeck_tensor = as_tensor([[S_xx, 0], [0, S_yy]])
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
(...)
problem = NonlinearProblem(weak_form, TempVolt, bcs=bcs, J = Jac)
(...)
n, converged = solver.solve(TempVolt)

# Current density.
J = VectorFunctionSpace(mesh, ("DG", deg-1))
Je = Function(J)
Je_flux_calculator = Expression(-sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp), J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)

# Curl of current density.
Jcurl = functionspace(mesh, ("DG", deg-2, (mesh.geometry.dim,)))
Jecurl = Function(Jcurl)
Jecurl_flux_calculator = Expression(curl(as_vector([J_vector[0],J_vector[1],0])), Jcurl.element.interpolation_points())
Jecurl.interpolate(Jecurl_flux_calculator)

```

with traceback:

```auto
Traceback (most recent call last):
  File "/root/shared/TE_simulation/Bridgman_Nye_example.py", line 434, in <module>
    post_proc.post_proc()
  File "/root/shared/TE_simulation/Bridgman_Nye_example.py", line 363, in post_proc
    J_curl = build_vector_field(curl(as_vector([J_vector[0],J_vector[1],0])), mesh, deg)
  File "/root/shared/TE_simulation/Bridgman_Nye_example.py", line 51, in build_vector_field
    flux_calculator = Expression(expression, Q.element.interpolation_points())
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 149, in __init__
    self._ufcx_expression, module, self._code = jit.ffcx_jit(comm, (e, _X),
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 56, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 210, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 247, in compile_expressions
    raise e
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 238, in compile_expressions
    impl = _compile_objects(decl, expressions, expr_names, module_name, p, cache_dir,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 260, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/compiler.py", line 97, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, options)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 95, in analyze_ufl_objects
    processed_expression = _analyze_expression(original_expression, options)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 117, in _analyze_expression
    expression = ufl.algorithms.apply_algebra_lowering.apply_algebra_lowering(expression)
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 148, in apply_algebra_lowering
    return map_integrand_dags(LowerCompoundAlgebra(), expr)
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 73, in map_integrand_dags
    return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 66, in map_integrands
    return function(integrand)
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 73, in <lambda>
    return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
  File "/usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py", line 34, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress,
  File "/usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py", line 100, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 142, in curl
    return as_vector((c(1, 2), c(2, 0), c(0, 1)))
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 135, in c
    return a[j].dx(i) - a[i].dx(j)
  File "/usr/local/lib/python3.10/dist-packages/ufl/exproperators.py", line 435, in _dx
    d = Grad(d)
  File "/usr/local/lib/python3.10/dist-packages/ufl/differentiation.py", line 252, in __new__
    dim = find_geometric_dimension(f)
  File "/usr/local/lib/python3.10/dist-packages/ufl/domain.py", line 396, in find_geometric_dimension
    raise ValueError("Cannot determine geometric dimension from expression.")
ValueError: Cannot determine geometric dimension from expression.

```

I also tried the Expression: `J_vector[1].dx(0)-J_vector[0].dx(1)` and `as_vector([0,0,J_vector[1].dx(0)-J_vector[0].dx(1)])`, without any success.

I am sorry for not posting my whole code (I can, but it’s hard to parse, I’d rather make a MWE).

And here’s my attempt to the MWE:

```auto
import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
                         dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
                         locate_dofs_topological, Expression, create_matrix, assemble_matrix)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
                 TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
                 inv, split, derivative, FacetNormal, div, curl, as_vector)

from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8    
T_cold = 100
    
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

def boundary_D_left(x):
    return np.isclose(x[0], 0)

def boundary_D_right(x):
    return np.isclose(x[0], 1)

# Define ME function space
deg = 3
el = FiniteElement("CG", mesh.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(mesh, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

# The Seebeck coefficient has to be an anisotropic tensor for the Bridgman heat to take place.
s_xx, s_yy = seebeck_components
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])

# Define the boundary conditions
#left_facets = facet_markers.find(inj_current_curve)
#right_facets = facet_markers.find(out_current_curve)
#left_dofs = locate_dofs_topological(ME.sub(1), mesh.topology.dim-1, left_facets)

#left_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim-1, left_facets)
#right_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim-1, right_facets)
left_dofs_temp = locate_dofs_geometrical(ME.sub(0), boundary_D_left)
right_dofs_temp = locate_dofs_geometrical(ME.sub(0), boundary_D_right)

bc_temp_left = dirichletbc(ScalarType(T_cold), left_dofs_temp, ME.sub(0))
bc_temp_right = dirichletbc(ScalarType(T_cold + ΔT), right_dofs_temp, ME.sub(0))
bcs = [bc_temp_left, bc_temp_right]

x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)
ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)

J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))

# Weak form.
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
Fourier_term = dot(-κ * grad(temp), grad(u)) * dx
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve) 

weak_form = F_T + F_V

# Solve the PDE.
Jac = derivative(weak_form,TempVolt,dTV)

print(f''' ------- Pre-processing --------
    Current density injected: {J_inj}
    Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))}
    Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))}
    This should correspond to the current injected: {assemble_scalar(form(J_out*ds(out_current_curve)))}
''')

# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-13
solver.atol = 1e-13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)    
n, converged = solver.solve(TempVolt)
assert(converged)
n = FacetNormal(mesh)
    
    
# Current density.
J = functionspace(mesh, ("DG", deg-1, (mesh.geometry.dim,)))
Je = Function(J)
Je_flux_calculator = Expression(J_vector, J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)
        
Jcurl = functionspace(mesh, ("DG", deg-1, (mesh.geometry.dim,)))
Jecurl = Function(Jcurl)
Jecurl_flux_calculator = Expression(curl(J_vector), Jcurl.element.interpolation_points())
Jecurl.interpolate(Jecurl_flux_calculator)
            
J_curl = build_vector_field(as_vector([0,0,J_vector[1].dx(0)-J_vector[0].dx(1)]), mesh, deg)
J_curl = build_vector_field(J_vector[1].dx(0)-J_vector[0].dx(1), mesh, deg)
J_curl = build_vector_field(curl(as_vector([J_vector[0],J_vector[1],0])), mesh, deg)

```

But `AttributeError: module 'dolfinx.mesh' has no attribute 'ufl_cell'` so I’ll have to debug myself.

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 20, 2024, 8:46pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/3 "2024-06-20T20:46:21Z")

</div>

Here is my new attempt to make a working MWE:

```auto
import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
                         dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
                         locate_dofs_topological, Expression, create_matrix, assemble_matrix)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
                 TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
                 inv, split, derivative, FacetNormal, div, curl, as_vector)

from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8    
T_cold = 100
    
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

def boundary_D_left(x):
    return np.isclose(x[0], 0)

def boundary_D_right(x):
    return np.isclose(x[0], 1)

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 1e-3, 1e-4
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])

left_dofs_temp = locate_dofs_geometrical(ME.sub(0), boundary_D_left)
right_dofs_temp = locate_dofs_geometrical(ME.sub(0), boundary_D_right)

bc_temp_left = dirichletbc(ScalarType(T_cold), left_dofs_temp, ME.sub(0))
bc_temp_right = dirichletbc(ScalarType(T_cold + ΔT), right_dofs_temp, ME.sub(0))
bcs = [bc_temp_left, bc_temp_right]

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain, subdomain_data=facet_markers)

J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve, domain=domain)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve, domain=domain)))

# Weak form.
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
Fourier_term = dot(-κ * grad(temp), grad(u)) * dx
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve) 

weak_form = F_T + F_V

Jac = derivative(weak_form,TempVolt,dTV)

# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-13
solver.atol = 1e-13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)    
n, converged = solver.solve(TempVolt)
assert(converged)
n = FacetNormal(domain)
    
    
# Current density.
J = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Je = Function(J)
Je_flux_calculator = Expression(J_vector, J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)
        
Jcurl = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Jecurl = Function(Jcurl)
Jecurl_flux_calculator = Expression(curl(J_vector), Jcurl.element.interpolation_points())
Jecurl.interpolate(Jecurl_flux_calculator)

```

The traceback:

```auto
Traceback (most recent call last):
  File "/root/shared/debug_discourse/curl_MWE.py", line 56, in <module>
    left_dofs_temp = locate_dofs_geometrical(ME.sub(0), boundary_D_left)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/bcs.py", line 51, in locate_dofs_geometrical
    return _cpp.fem.locate_dofs_geometrical(V._cpp_object, marker) # type: ignore
RuntimeError: Cannot tabulate coordinates for a FunctionSpace that is a subspace.

```

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [June 21, 2024, 4:40am UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/4 "2024-06-21T04:40:54Z")

</div>

> [@raboynics](#):
>
> ```auto
> left_dofs_temp = locate_dofs_geometrical(ME.sub(0), boundary_D_left)
> 
> ```

Try to follow one of the demos that uses `mixed_element`, for instance [dolfinx/python/demo/demo\_stokes.py at v0.8.0 · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/v0.8.0/python/demo/demo_stokes.py#L474-L479) . (the fact that `locate_dofs_topological` is used there instead of `locate_dofs_geometrical` shouldn’t make any difference)

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 21, 2024, 8:22pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/5 "2024-06-21T20:22:26Z")

</div>

Grazie Francesco, I am able to make some progress thanks to you.

I am a bit stuck. My latest attempt:

```auto
import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
                         dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
                         locate_dofs_topological, Expression, create_matrix, assemble_matrix,
                         )
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
                 TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
                 inv, split, derivative, FacetNormal, div, curl, as_vector)
from dolfinx.mesh import locate_entities_boundary
from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8    
T_cold = 100
 
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
def boundary_D_left(x):
    return np.isclose(x[0], 0)

def boundary_D_right(x):
    return np.isclose(x[0], 1)

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 1e-3, 1e-4
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])

T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
facets = locate_entities_boundary(domain, 1, boundary_D_right)
dofs = locate_dofs_topological((T0, Q), 1, facets)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T_cold = Function(Q)
facets = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T0, Q), 1, facets)
bc_temp_right = dirichletbc(T_cold, dofs, T0)

bcs = [bc_temp_left, bc_temp_right]

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain)#, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain)#, subdomain_data=facet_markers)

the_current = 1.0
J_inj = the_current / assemble_scalar(form(1 * ds))
J_out = the_current / assemble_scalar(form(1 * ds))

```

1st question: I do not understand in the example where is the value of “noslip” specified. I had T\_cold and T\_hot as floats before (or Constants), but now they are “Function(Q)” without any value. How do I specify their values? I hope not by interpolating a constant into Q?

2nd question: How do I specify ds(1) and ds(2) for example for the 2 surfaces or curves where I apply the Dirichlet b.c.? I would need to mark the curves? This is not done in the example? I am confused on this.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [June 22, 2024, 4:51am UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/6 "2024-06-22T04:51:44Z")

</div>

> [@raboynics](#):
>
> 1st question: I do not understand in the example where is the value of “noslip” specified. I had T\_cold and T\_hot as floats before (or Constants), but now they are “Function(Q)” without any value. How do I specify their values? I hope not by interpolating a constant into Q?

A `dolfinx.fem.Function` has zero values by default, which is indeed the value for no split boundary conditions. Interpolating the constant value into a CG space is simply setting all entries of the `dolfinx.fem.Function.x.array` to that value. If someone else reads your code, you may want to be more expressive and interpolate.

> [@raboynics](#):
>
> 2nd question: How do I specify ds(1) and ds(2) for example for the 2 surfaces or curves where I apply the Dirichlet b.c.? I would need to mark the curves? This is not done in the example? I am confused on this.

Since you are already locating the entities, you may want to use a `meshtags` defined on those entities as in [dolfinx/python/test/unit/fem/test\_assemble\_domains.py at v0.8.0 · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/v0.8.0/python/test/unit/fem/test_assemble_domains.py#L270-L274)  
In principle you could even avoid the `meshtags` object if adopting the strategy in [dolfinx/python/test/unit/fem/test\_assemble\_domains.py at v0.8.0 · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/v0.8.0/python/test/unit/fem/test_assemble_domains.py#L307-L318) . I’ll leave it to you to decide what you prefer.

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 22, 2024, 9:29pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/7 "2024-06-22T21:29:27Z")

</div>

Thanks once more Francesco.  
I could proceed further with my MWE, but my code fails to solve the problem (Newton method fails to converge, and in fact diverges at the very first step).  
I tried to simplify the equation (made S tensor equal to 0, removed Fourier and Bridgman terms in the weak form), instead of -nan I sometimes got “inf” according to the logger.  
The traceback is:

```auto
area of material: 1.0
1.0000000000000002
1.0000000000000002
 ------- Pre-processing --------
    Current density injected: 1.0000000000000002
    Length of the side where current is injected: 0.9999999999999999
    Length of the side where current leaves the wire: 0.9999999999999999
    This should correspond to the current injected: 0.9999999999999999

2024-06-22 21:22:24.037 ( 0.740s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 0: r (abs) = -nan (tol = 1e-06) r (rel) = -nan(tol = 1e-06)
2024-06-22 21:22:24.039 ( 0.743s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.045 ( 0.749s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.005737, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.046 ( 0.749s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 1: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.048 ( 0.751s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.051 ( 0.754s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.003069, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.052 ( 0.755s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 2: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.054 ( 0.757s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.054 ( 0.757s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000033, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.055 ( 0.758s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 3: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.057 ( 0.760s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.057 ( 0.760s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000047, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.058 ( 0.761s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 4: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.060 ( 0.763s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.060 ( 0.763s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000021, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.060 ( 0.764s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 5: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.062 ( 0.766s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.062 ( 0.766s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000018, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.063 ( 0.766s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 6: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.065 ( 0.768s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.065 ( 0.768s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000018, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.066 ( 0.769s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 7: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.068 ( 0.771s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.068 ( 0.771s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000021, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.068 ( 0.772s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 8: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.070 ( 0.774s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.070 ( 0.774s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000018, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.071 ( 0.774s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 9: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.073 ( 0.777s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.073 ( 0.777s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000029, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.074 ( 0.777s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 10: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.076 ( 0.779s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.076 ( 0.780s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000038, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.077 ( 0.780s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 11: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.079 ( 0.782s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.079 ( 0.782s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000018, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.079 ( 0.783s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 12: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.081 ( 0.785s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.081 ( 0.785s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000017, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.082 ( 0.785s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 13: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.084 ( 0.787s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.084 ( 0.788s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000018, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.085 ( 0.788s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 14: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.087 ( 0.790s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.087 ( 0.790s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000017, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.087 ( 0.791s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 15: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.089 ( 0.793s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.089 ( 0.793s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000039, 0.000000, 0.000000 (PETSc Krylov solver)
(...)
2024-06-22 21:22:24.317 ( 1.020s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000021, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.317 ( 1.021s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 99: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.319 ( 1.023s) [main] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.319 ( 1.023s) [main] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000022, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.320 ( 1.023s) [main] NewtonSolver.cpp:36 INFO| Newton iteration 100: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
Traceback (most recent call last):
  File "/root/shared/debug_discourse/curl_MWE.py", line 153, in <module>
    n, converged = solver.solve(TempVolt)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/nls/petsc.py", line 46, in solve
    n, converged = super().solve(u.vector)
RuntimeError: Newton solver did not converge because maximum number of iterations reached
2024-06-22 21:22:24.366 ( 1.069s) [main] loguru.cpp:526 INFO| atexit

```

My full code:

```auto
import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
                         dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
                         locate_dofs_topological, Expression, create_matrix, assemble_matrix,
                         )
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
                 TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
                 inv, split, derivative, FacetNormal, div, curl, as_vector)
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8    
    
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
x = SpatialCoordinate(domain)

def boundary_D_left(x):
    return np.isclose(x[0], 0)

def boundary_D_right(x):
    return np.isclose(x[0], 1)

def T_cold_func(x): 
    return eval(str(np.asarray([100.0])))
def T_hot_func(x): 
    return eval(str(np.asarray([200.0])))

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 0,0#1e-3, 1e-4
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])

T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
T_hot.interpolate(T_hot_func)
facets_hot = locate_entities_boundary(domain, 1, boundary_D_right)
dofs = locate_dofs_topological((T0, Q), 1, facets_hot)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T_cold = Function(Q)
T_cold.interpolate(T_cold_func)
facets_cold = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T0, Q), 1, facets_cold)
bc_temp_right = dirichletbc(T_cold, dofs, T0)

bcs = [bc_temp_left, bc_temp_right]

# Mesh tags to apply Neumann b.c.
tdim = domain.topology.dim
domain.topology.create_entities(tdim - 1)
facet_map = domain.topology.index_map(tdim - 1)
num_facets = facet_map.size_local + facet_map.num_ghosts
facet_indices = np.arange(0, num_facets)
facet_values = np.zeros_like(facet_indices, dtype=np.intc)

# T_cold surface markers.
inj_current_curve = 12
out_current_curve = 13
facet_values[facets_cold] = inj_current_curve
facet_values[facets_hot] = out_current_curve
mt_facets = meshtags(domain, tdim - 1, facet_indices, facet_values)

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain)#, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain, subdomain_data=mt_facets)

the_current = 1.0
J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve)))

print(f'area of material: {assemble_scalar(form(1 * dx))}')
print(J_inj)
print(J_out)
# Weak form.
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
Fourier_term = dot(-κ * grad(temp), grad(u)) * dx
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve) 

weak_form = F_T + F_V

# Solve the PDE.
Jac = derivative(weak_form,TempVolt,dTV)

print(f''' ------- Pre-processing --------
    Current density injected: {J_inj}
    Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=domain)))}
    Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=domain)))}
    This should correspond to the current injected: {assemble_scalar(form(J_out*ds(out_current_curve)))}
''')

# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-6#13
solver.atol = 1e-6#13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)    
n, converged = solver.solve(TempVolt)
assert(converged)
n = FacetNormal(domain)
    
    
# Current density.
J = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Je = Function(J)
Je_flux_calculator = Expression(J_vector, J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)
        
Jcurl = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Jecurl = Function(Jcurl)
Jecurl_flux_calculator = Expression(curl(J_vector), Jcurl.element.interpolation_points())
Jecurl.interpolate(Jecurl_flux_calculator)
            
J_curl = build_vector_field(as_vector([0,0,J_vector[1].dx(0)-J_vector[0].dx(1)]), domain, deg)
J_curl = build_vector_field(J_vector[1].dx(0)-J_vector[0].dx(1), domain, deg)
J_curl = build_vector_field(curl(as_vector([J_vector[0],J_vector[1],0])), domain, deg)

```

I’m very close to have it working, which would be awesome… But not sure how to make it converge now.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [June 23, 2024, 5:42am UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/8 "2024-06-23T05:42:44Z")

</div>

An initial residual with value `nan` or `inf` typically means that the initial guess you provided to the Newton solver, i.e. the value of `TempVolt`, is incorrect. Unless you interpolate some expression in `TempVolt` is zero, which may lead to `inf`/`nan` as soon as your residual has expressions like `1/TempVolt` (or similar algebraic operations that do not make sense when the solution is zero).

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 23, 2024, 8:38pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/9 "2024-06-23T20:38:12Z")

</div>

I see…  
I’m having a hard time to do that, i.e. interpolate an intial guess for temp and volt.  
Here is my attempt (exact same whole code as above, with this addition):

```auto
def initial_guess_temp():
    return np.asarray([100.0])
    
def initial_guess_volt():
    return np.asarray([0.0])
TempVolt.sub(0).interpolate(initial_guess_temp, facet_indices)
TempVolt.sub(1).interpolate(initial_guess_volt, facet_indices)
# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)

```

This yields:

```auto
area of material: 1.0
0.0
0.0
 ------- Pre-processing --------
    Current density injected: 0.0
    Length of the side where current is injected: 0.9999999999999999
    Length of the side where current leaves the wire: 0.9999999999999999
    This should correspond to the current injected: 0.0

Traceback (most recent call last):
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 397, in interpolate
    _interpolate(u, cells)
  File "/usr/lib/python3.10/functools.py", line 889, in wrapper
    return dispatch(args[0]. __class__ )(*args, **kw)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-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 0x7da93e7f9770>, <function initial_guess_temp at 0x7da93e4ebf40>, 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, 64,
        65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,
        78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
        91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
       104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
       117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
       130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
       143]), ((), (), (), ())

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/root/shared/debug_discourse/curl_MWE.py", line 133, in <module>
    TempVolt.sub(0).interpolate(initial_guess_temp, facet_indices)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 401, in interpolate
    x = _cpp.fem.interpolation_coords(self._V.element, self._V.mesh.geometry, cells)
TypeError: interpolation_coords(): incompatible function arguments. The following argument types are supported:
    1. (element: dolfinx.cpp.fem.FiniteElement_float32, V: dolfinx.cpp.mesh.Geometry_float32, cells: numpy.ndarray[numpy.int32]) -> numpy.ndarray[numpy.float32]
    2. (element: dolfinx.cpp.fem.FiniteElement_float64, V: dolfinx.cpp.mesh.Geometry_float64, cells: numpy.ndarray[numpy.int32]) -> numpy.ndarray[numpy.float64]

Invoked with: <dolfinx.cpp.fem.FiniteElement_float64 object at 0x7da93e2ed3b0>, <dolfinx.cpp.mesh.Geometry_float64 object at 0x7da93e2eea70>, 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, 64,
        65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,
        78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
        91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
       104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
       117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
       130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
       143])

```

I also have a hard time to find ME examples with interpolation.

---

<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:** [June 23, 2024, 8:59pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/10 "2024-06-23T20:59:28Z")

</div>

> [@raboynics](#):
>
> ```auto
> TempVolt.sub(0).interpolate(initial_guess_temp, facet_indices)
> TempVolt.sub(1).interpolate(initial_guess_volt, facet_indices)
> 
> ```

You cannot interpolate on facets only.

What i would do is using `compute_incident_entities` to find all cells connected to `facet_indices` and send this into interpolate.

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 24, 2024, 9:37pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/11 "2024-06-24T21:37:11Z")

</div>

I see, thanks a lot for your help.  
What is wrong with my new attempt though?

```auto
def initial_guess_temp():
    return np.asarray([100.0])
    
def initial_guess_volt():
    return np.asarray([0.0])

the_cells = compute_incident_entities(domain.topology, facet_indices, tdim, tdim-1)

TempVolt.sub(0).interpolate(initial_guess_temp, the_cells)
TempVolt.sub(1).interpolate(initial_guess_volt, the_cells)

```

The output is:

```auto
area of material: 1.0
0.0
0.0
 ------- Pre-processing --------
    Current density injected: 0.0
    Length of the side where current is injected: 0.9999999999999999
    Length of the side where current leaves the wire: 0.9999999999999999
    This should correspond to the current injected: 0.0

Traceback (most recent call last):
  File "/root/shared/debug_discourse/curl_MWE.py", line 134, in <module>
    the_cells = compute_incident_entities(domain.topology, facet_indices, tdim, tdim-1)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/mesh.py", line 182, in compute_incident_entities
    return _cpp.mesh.compute_incident_entities(topology, entities, d0, d1)
TypeError: compute_incident_entities(): incompatible function arguments. The following argument types are supported:
    1. (mesh: dolfinx.cpp.mesh.Topology, entities: numpy.ndarray[numpy.int32], d0: int, d1: int) -> numpy.ndarray[numpy.int32]

Invoked with: <dolfinx.cpp.mesh.Topology object at 0x7c21e66ac220>, 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, 64,
        65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,
        78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
        91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
       104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
       117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
       130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
       143]), 2, 1

```

---

<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:** [June 25, 2024, 5:48am UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/12 "2024-06-25T05:48:03Z")

</div>

Please note that it is hard to give you guidance when you provide small code snippets that are not reproducible.  
Say, for the problem above, you should start by testing the method on a unit square.

```python
from mpi4py import MPI
import dolfinx
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
v = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
u = dolfinx.fem.Function(v)

def top_facets(x):
    return np.isclose(x[1], 1.0)

tdim = mesh.topology.dim
fdim = tdim -1

mesh.topology.create_connectivity(fdim, tdim)
facets = dolfinx.mesh.locate_entities(mesh, fdim, top_facets)

cells = dolfinx.mesh.compute_incident_entities(mesh.topology, facets, fdim, tdim)

u.interpolate(lambda x: np.sin(x[0]), cells=cells)
u.x.scatter_forward()

with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u], engine="BP4") as bp:
    bp.write(0.0)

```

returning

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/0/0d46bbf4ab9fd69993099c35f3e596e4aa55f275.png)

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 25, 2024, 3:16pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/13 "2024-06-25T15:16:44Z")

</div>

> [@dokken](#):
>
> ```auto
> tdim = mesh.topology.dim
> fdim = tdim -1
> 
> mesh.topology.create_connectivity(fdim, tdim)
> facets = dolfinx.mesh.locate_entities(mesh, fdim, top_facets)
> 
> cells = dolfinx.mesh.compute_incident_entities(mesh.topology, facets, fdim, tdim)
> 
> ```

I’m sorry dokken, even with your code, I am having a hard time. I do not want to do exactly what you’ve done, I believe. I want to set “temp” and “volt” to constant values over the whole mesh, not on particular submeshes.  
I am not sure, maybe I am missing the `domain.topology.create_connectivity(tdim-1, tdim)` part (I tried to add it but this didn’t fix my error), but I am not sure because I do not see it in the code linked by Francesco. Also, you seem to have altered the order of tdim and fdim in the arguments of `compute_incident_entities`. Which order would be correct?

Here is my full code so far:

```auto
import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
                         dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
                         locate_dofs_topological, Expression, create_matrix, assemble_matrix,
                         )
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
                 TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
                 inv, split, derivative, FacetNormal, div, curl, as_vector)
from dolfinx.mesh import locate_entities_boundary, meshtags, compute_incident_entities
from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8    
    
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
x = SpatialCoordinate(domain)

def boundary_D_left(x):
    return np.isclose(x[0], 0)

def boundary_D_right(x):
    return np.isclose(x[0], 1)

def T_cold_func(x): 
    return eval(str(np.asarray([100.0])))
def T_hot_func(x): 
    return eval(str(np.asarray([100.0])))

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 0,0
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])
the_current = 0.0

T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
T_hot.interpolate(T_hot_func)
facets_hot = locate_entities_boundary(domain, 1, boundary_D_right)
dofs = locate_dofs_topological((T0, Q), 1, facets_hot)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T_cold = Function(Q)
T_cold.interpolate(T_cold_func)
facets_cold = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T0, Q), 1, facets_cold)
bc_temp_right = dirichletbc(T_cold, dofs, T0)

bcs = [bc_temp_left, bc_temp_right]

# Mesh tags to apply Neumann b.c.
tdim = domain.topology.dim
domain.topology.create_entities(tdim - 1)
facet_map = domain.topology.index_map(tdim - 1)
num_facets = facet_map.size_local + facet_map.num_ghosts
facet_indices = np.arange(0, num_facets)
facet_values = np.zeros_like(facet_indices, dtype=np.intc)

# T_cold surface markers.
inj_current_curve = 12
out_current_curve = 13
facet_values[facets_cold] = inj_current_curve
facet_values[facets_hot] = out_current_curve
mt_facets = meshtags(domain, tdim - 1, facet_indices, facet_values)

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain)#, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain, subdomain_data=mt_facets)

J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve)))

print(f'area of material: {assemble_scalar(form(1 * dx))}')
print(J_inj)
print(J_out)
# Weak form.
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
Fourier_term = dot(-κ * grad(temp), grad(u)) * dx
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve) 

weak_form = F_T + F_V

# Solve the PDE.
Jac = derivative(weak_form,TempVolt,dTV)

print(f''' ------- Pre-processing --------
    Current density injected: {J_inj}
    Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=domain)))}
    Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=domain)))}
    This should correspond to the current injected: {assemble_scalar(form(J_out*ds(out_current_curve)))}
''')

#print(facet_indices)
# Create initial guess.
def initial_guess_temp():
    return np.asarray([100.0])
    
def initial_guess_volt():
    return np.asarray([0.0])
    
tdim = mesh.topology.dim
fdim = tdim -1 
    
domain.topology.create_connectivity(tdim-1, tdim)
the_cells = compute_incident_entities(domain.topology, facet_indices, tdim-1, tdim)

TempVolt.sub(0).interpolate(initial_guess_temp, the_cells)
TempVolt.sub(1).interpolate(initial_guess_volt, the_cells)
# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"#"incremental"#"residual"
solver.rtol = 1e-6#13
solver.atol = 1e-6#13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)    
n, converged = solver.solve(TempVolt)
assert(converged)
n = FacetNormal(domain)

```

---

<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:** [June 25, 2024, 3:21pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/14 "2024-06-25T15:21:22Z")

</div>

> [@raboynics](#):
>
> I’m sorry dokken, even with your code, I am having a hard time. I do not want to do exactly what you’ve done, I believe. I want to set “temp” and “volt” to constant values over the whole mesh, not on particular submeshes.

But why dont you the just do interpolation over the whole mesh.

I suggested compute\_incident\_entities because you had

> [@raboynics](#):
>
> ```auto
> def initial_guess_temp():
> return np.asarray([100.0])
>     
> def initial_guess_volt():
> return np.asarray([0.0])
> TempVolt.sub(0).interpolate(initial_guess_temp, facet_indices)
> TempVolt.sub(1).interpolate(initial_guess_volt, facet_indices)
> # Solve the PDE.
> problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
> 
> ```

in your code, which is simply not correct.

What is wrong with doing  
`TempVolt.sub(0).interpolate(initial_guess_temp)`?

> [@raboynics](#):
>
> Also, you seem to have altered the order of tdim and fdim in the arguments of `compute_incident_entities`. Which order would be correct?

The order I have in my code is correct. First dimension is the dimension of the input entities, second dimension is output dimension.  
See for instance: [What is the difference between compute\_incident\_entities and entities\_to\_geometry?](https://fenicsproject.discourse.group/t/what-is-the-difference-between-compute-incident-entities-and-entities-to-geometry/14867)

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 25, 2024, 3:37pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/15 "2024-06-25T15:37:59Z")

</div>

> [@dokken](#):
>
> But why dont you the just do interpolation over the whole mesh.
> 
> I suggested compute\_incident\_entities because you had

I guess that’s what I tried and what I thought I was doing. I did not realize I wasn’t doing this.

> [@dokken](#):
>
> What is wrong with doing  
> `TempVolt.sub(0).interpolate(initial_guess_temp)`?

Answer: `TypeError: initial_guess_temp() takes 0 positional arguments but 1 was given`. This was the very first thing I tried.  
The full traceback is:

```auto
Traceback (most recent call last):
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 397, in interpolate
    _interpolate(u, cells)
  File "/usr/lib/python3.10/functools.py", line 889, in wrapper
    return dispatch(args[0]. __class__ )(*args, **kw)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-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 0x7b0a0c6ce770>, <function initial_guess_temp at 0x7b0a0c21fa30>, 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 "/root/shared/debug_discourse/curl_MWE.py", line 127, in <module>
    TempVolt.sub(0).interpolate(initial_guess_temp)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 402, in interpolate
    self._cpp_object.interpolate(np.asarray(u(x), dtype=self.dtype), cells)
TypeError: initial_guess_temp() takes 0 positional arguments but 1 was given

```

> [@dokken](#):
>
> The order I have in my code is correct. First dimension is the dimension of the input entities, second dimension is output dimension.

Great. Note that this differs from your Github issue (I’m sure you’re aware… but just in case): [[BUG]: Packing of codim-0 form not working for mixture of measures · Issue #3256 · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/issues/3256).

---

<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:** [June 25, 2024, 3:43pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/16 "2024-06-25T15:43:47Z")

</div>

> [@raboynics](#):
>
> Great. Note that this differs from your Github issue (I’m sure you’re aware… but just in case): [[BUG]: Packing of codim-0 form not working for mixture of measures · Issue #3256 · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/issues/3256).

Thats because there i am mapping cells to facets, while you are mapping facets to cells.

Ref your function:

> [@raboynics](#):
>
> ```auto
> def initial_guess_temp():
> return np.asarray([100.0])
>     
> def initial_guess_volt():
> return np.asarray([0.0])
> 
> ```

Change them to

```auto
def initial_guess_temp(x):
    return np.full(x.shape[1], 100.)
    
def initial_guess_volt(x):
    return np.zeros(x.shape[1])

```

The interpolate function takes in an x that is (3, num\_points) and returns an array of length num\_points mapping each entry x[:,i] to a corresponding output value

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 25, 2024, 3:56pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/17 "2024-06-25T15:56:57Z")

</div>

Thanks a lot dokken. I had found an older code of mine and the syntax was a bit different, I couldn’t make it work, but your code did the trick.

I still get the same problem of instant nan values with Newton method. I am guessing something else is wrong in my code.

I might try to impose a Dirichlet b.c. for “volt” although I do not have it in my original (non MWE) code.

I do have a question regarding the spaces of the boundary condition functions:  
Should I create different spaces “T0” and “T1” as well as “Q” and “Q1” like so, one for each of the Dirichlet b.c. for the “temp” variable, like so?:

```auto
# Define the boundary conditions
T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
T_hot.interpolate(T_hot_func)
facets_hot = locate_entities_boundary(domain, 1, boundary_D_right)
dofs = locate_dofs_topological((T0, Q), 1, facets_hot)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T1 = ME.sub(0)
Q1, _ = T1.collapse()
T_cold = Function(Q1)
T_cold.interpolate(T_cold_func)
facets_cold = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T1, Q1), 1, facets_cold)
bc_temp_right = dirichletbc(T_cold, dofs, T1)

```

Or is the code I posted previously OK regarding this part?

Either way, this does not fix Newton method failing to converge.

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 25, 2024, 4:08pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/18 "2024-06-25T16:08:10Z")

</div>

Here is my latest code with the implementation of a Dirichlet b.c. for volt, Newton method fails to converge:

```auto
import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
                         dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
                         locate_dofs_topological, Expression, create_matrix, assemble_matrix,
                         )
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
                 TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
                 inv, split, derivative, FacetNormal, div, curl, as_vector)
from dolfinx.mesh import locate_entities_boundary, meshtags, compute_incident_entities
from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8    
    
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
x = SpatialCoordinate(domain)

def boundary_D_left(x):
    return np.isclose(x[0], 0)

def boundary_D_right(x):
    return np.isclose(x[0], 1)

def T_cold_func(x): 
    return eval(str(np.asarray([100.0])))
def T_hot_func(x): 
    return eval(str(np.asarray([100.0])))
    
def V_left_func(x):
    return eval(str(np.asarray([0.0])))

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 0,0#1e-3, 1e-4
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])
the_current = 0.0

# Define the boundary conditions
T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
T_hot.interpolate(T_hot_func)
facets_hot = locate_entities_boundary(domain, 1, boundary_D_right)
dofs = locate_dofs_topological((T0, Q), 1, facets_hot)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T1 = ME.sub(0)
Q1, _ = T1.collapse()
T_cold = Function(Q1)
T_cold.interpolate(T_cold_func)
facets_cold = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T1, Q1), 1, facets_cold)
bc_temp_right = dirichletbc(T_cold, dofs, T1)

V0 = ME.sub(1)
JV, _ = V0.collapse()
V_hot = Function(JV)
V_hot.interpolate(V_left_func)
bc_volt_left = dirichletbc(V_hot, dofs, V0)

bcs = [bc_temp_left, bc_temp_right, bc_volt_left]

# Mesh tags to apply Neumann b.c.
tdim = domain.topology.dim
fdim = tdim -1 
domain.topology.create_entities(fdim)
facet_map = domain.topology.index_map(fdim)
num_facets = facet_map.size_local + facet_map.num_ghosts
facet_indices = np.arange(0, num_facets)
facet_values = np.zeros_like(facet_indices, dtype=np.intc)

# T_cold surface markers.
inj_current_curve = 12
out_current_curve = 13
facet_values[facets_cold] = inj_current_curve
facet_values[facets_hot] = out_current_curve
mt_facets = meshtags(domain, fdim, facet_indices, facet_values)

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain)#, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain, subdomain_data=mt_facets)

J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve)))

print(f'area of material: {assemble_scalar(form(1 * dx))}')
print(J_inj)
print(J_out)
# Weak form.
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
Fourier_term = dot(-κ * grad(temp), grad(u)) * dx
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve) 

weak_form = F_T + F_V

# Solve the PDE.
Jac = derivative(weak_form,TempVolt,dTV)

print(f''' ------- Pre-processing --------
    Current density injected: {J_inj}
    Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=domain)))}
    Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=domain)))}
    This should correspond to the current injected: {assemble_scalar(form(J_out*ds(out_current_curve)))}
''')

def initial_guess_temp(x):
    return np.full(x.shape[1], 100.0)
    
def initial_guess_volt(x):
    return np.zeros(x.shape[1])
         
TempVolt.sub(0).interpolate(initial_guess_temp)
TempVolt.sub(1).interpolate(initial_guess_volt)

# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"#"incremental"#"residual"
solver.rtol = 1e-6#13
solver.atol = 1e-6#13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)    
n, converged = solver.solve(TempVolt)
assert(converged)   

```

---

<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:** [June 25, 2024, 5:45pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/19 "2024-06-25T17:45:34Z")

</div>

Consider the following:

```python
import dolfinx.fem
vec = dolfinx.fem.assemble_vector(dolfinx.fem.form(weak_form))
dolfinx.fem.apply_lifting(vec.array, [dolfinx.fem.form(Jac)], [bcs])

```

This would be part of your RHS of the newton problem (with some extra stuff).  
Inspecting `vec.array` yields:

```python
array([-7.54294650e+307, -2.08755243e+307, -2.02282709e+307, ...,
        0.00000000e+000, 0.00000000e+000, 0.00000000e+000])

```

indicating that something is seriously wrong with either your jacobian or your bcs.  
Further investigating, you can compute:

```python
A = dolfinx.fem.assemble_matrix(dolfinx.fem.form(Jac))
print(A.to_scipy().min(), A.to_scipy().max())

```

and see that the values are bounded.  
Lets next consider your DirichletBC.

`V_hot.x.array` contains `nan`’s and infinite values, meaning that it is wrong.  
I’ve now given you a full recipe on how I address convergence issues, which you should take into account for further questions.

I have no clue why you are creating interpolation functions like this

> [@raboynics](#):
>
> ```auto
> def T_cold_func(x): 
> return eval(str(np.asarray([100.0])))
> def T_hot_func(x): 
> return eval(str(np.asarray([100.0])))
>     
> def V_left_func(x):
> return eval(str(np.asarray([0.0])))
> 
> ```

as I have explicitly stated how to create them in my last post on this thread.

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 25, 2024, 8:30pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/20 "2024-06-25T20:30:22Z")

</div>

Thank you very much dokken, I was sloppy, and unaware of these techniques to debug. This will help me in the future.  
I am very happy in that I reach the same point as with my original code (not the MWE).  
The code solves for temp and volt and computes “Je”, the current density of electrodynamics. This vector field, at first glance, looks correct to me (it has vortices in this case).  
I wish to compute its curl, and this is where I get stuck. Namely, the very last line of my script yields an error.  
Here’s the full code, the traceback and a picture of the current density vector field whose curl I am looking for. Note that I removed the Dirichlet b.c. on volt, the code is able to converge with pure Neumann b.c. and no null space specified.

```auto
import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
                         dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
                         locate_dofs_topological, Expression, create_matrix, assemble_matrix,
                         assemble_vector, apply_lifting)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
                 TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
                 inv, split, derivative, FacetNormal, div, curl, as_vector)
from dolfinx.mesh import locate_entities_boundary, meshtags, compute_incident_entities
from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8    
    
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
x = SpatialCoordinate(domain)

def boundary_D_left(x):
    return np.isclose(x[0], 0)

def boundary_D_bottom(x):
    return np.isclose(x[1], 0)

def T_cold_func(x): 
    return np.full(x.shape[1], 100.0)
def T_hot_func(x): 
    return np.full(x.shape[1], 200.0)
    
# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 1e-3, 1e-4
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])
the_current = 0.0

# Define the boundary conditions
T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
T_hot.interpolate(T_hot_func)
facets_hot = locate_entities_boundary(domain, 1, boundary_D_bottom)
dofs = locate_dofs_topological((T0, Q), 1, facets_hot)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T1 = ME.sub(0)
Q1, _ = T1.collapse()
T_cold = Function(Q1)
T_cold.interpolate(T_cold_func)
facets_cold = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T1, Q1), 1, facets_cold)
bc_temp_right = dirichletbc(T_cold, dofs, T1)
bcs = [bc_temp_left, bc_temp_right]#, bc_volt_left]

tdim = domain.topology.dim
fdim = tdim -1 
domain.topology.create_entities(fdim)
facet_map = domain.topology.index_map(fdim)
num_facets = facet_map.size_local + facet_map.num_ghosts
facet_indices = np.arange(0, num_facets)
facet_values = np.zeros_like(facet_indices, dtype=np.intc)

inj_current_curve = 12
out_current_curve = 13
facet_values[facets_cold] = inj_current_curve
facet_values[facets_hot] = out_current_curve
mt_facets = meshtags(domain, fdim, facet_indices, facet_values)

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain)#, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain, subdomain_data=mt_facets)

J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve)))

J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
Fourier_term = dot(-κ * grad(temp), grad(u)) * dx
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve) 

weak_form = F_T + F_V

Jac = derivative(weak_form,TempVolt,dTV)

def initial_guess_temp(x):
    return np.full(x.shape[1], 100.0)
    
def initial_guess_volt(x):
    return np.zeros(x.shape[1])
         
TempVolt.sub(0).interpolate(initial_guess_temp)
TempVolt.sub(1).interpolate(initial_guess_volt)

# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"#"incremental"#"residual"
solver.rtol = 1e-8#13
solver.atol = 1e-8#13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)    
n, converged = solver.solve(TempVolt)
assert(converged)
n = FacetNormal(domain)
    
# Current density.
J = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Je = Function(J)
Je_flux_calculator = Expression(J_vector, J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)

with VTXWriter(MPI.COMM_WORLD, "results/Je_debug.bp", Je, engine="BP4") as vtx:
    vtx.write(0.0)
    
Jcurl = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Jecurl = Function(Jcurl)
Jecurl_flux_calculator = Expression(curl(J_vector), Jcurl.element.interpolation_points())
Jecurl.interpolate(Jecurl_flux_calculator)
            
# Also tried:
#Jecurl_flux_calculator = Expression(as_vector([0,0,J_vector[1].dx(0)-J_vector[0].dx(1)]), Jcurl.element.interpolation_points())
#Jecurl_flux_calculator = Expression(J_vector[1].dx(0)-J_vector[0].dx(1), Jcurl.element.interpolation_points())
#Jecurl_flux_calculator = Expression(curl(as_vector([J_vector[0],J_vector[1],0])), Jcurl.element.interpolation_points())

```

```auto
Traceback (most recent call last):
  File "/root/shared/debug_discourse/curl_MWE.py", line 181, in <module>
    Jecurl.interpolate(Jecurl_flux_calculator)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 397, in interpolate
    _interpolate(u, cells)
  File "/usr/lib/python3.10/functools.py", line 889, in wrapper
    return dispatch(args[0]. __class__ )(*args, **kw)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 388, in _
    self._cpp_object.interpolate(expr._cpp_object, cells)
RuntimeError: Function value size not equal to Expression value size

```

The last commented lines all yield the same error except the very last, whose traceback is:

```auto
Traceback (most recent call last):
  File "/root/shared/debug_discourse/curl_MWE.py", line 183, in <module>
    Jecurl_flux_calculator = Expression(curl(as_vector([J_vector[0],J_vector[1],0])), Jcurl.element.interpolation_points())
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 149, in __init__
    self._ufcx_expression, module, self._code = jit.ffcx_jit(comm, (e, _X),
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 56, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 210, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 247, in compile_expressions
    raise e
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 238, in compile_expressions
    impl = _compile_objects(decl, expressions, expr_names, module_name, p, cache_dir,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 260, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/compiler.py", line 97, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, options)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 95, in analyze_ufl_objects
    processed_expression = _analyze_expression(original_expression, options)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 117, in _analyze_expression
    expression = ufl.algorithms.apply_algebra_lowering.apply_algebra_lowering(expression)
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 148, in apply_algebra_lowering
    return map_integrand_dags(LowerCompoundAlgebra(), expr)
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 73, in map_integrand_dags
    return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 66, in map_integrands
    return function(integrand)
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 73, in <lambda>
    return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
  File "/usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py", line 34, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress,
  File "/usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py", line 100, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 142, in curl
    return as_vector((c(1, 2), c(2, 0), c(0, 1)))
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 135, in c
    return a[j].dx(i) - a[i].dx(j)
  File "/usr/local/lib/python3.10/dist-packages/ufl/exproperators.py", line 435, in _dx
    d = Grad(d)
  File "/usr/local/lib/python3.10/dist-packages/ufl/differentiation.py", line 252, in __new__
    dim = find_geometric_dimension(f)
  File "/usr/local/lib/python3.10/dist-packages/ufl/domain.py", line 396, in find_geometric_dimension
    raise ValueError("Cannot determine geometric dimension from expression.")
ValueError: Cannot determine geometric dimension from expression.

```

 ![Je_debug](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/7/734c5c66c31278a5dffe9b9f22aac440d6204146.jpeg)

---

<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:** [June 25, 2024, 8:37pm UTC](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979/21 "2024-06-25T20:37:27Z")

</div>

You are using a 2D mesh.  
The ufl curl operator is only defined in 3D, as written in the documentation:  
[https://docs.fenicsproject.org/ufl/2024.1.0/manual/form\_language.html#curl-and-rot](https://docs.fenicsproject.org/ufl/2024.1.0/manual/form_language.html#curl-and-rot)

See for instance

> <https://github.com/FEniCS/dolfinx/blob/9f33dabc270b8e3c31e815515ecac5eb2a5414b0/python/demo/demo_pml.py#L378>

or [Vector curl in 2D - #2 by nate](https://fenicsproject.discourse.group/t/vector-curl-in-2d/5541/2)  
for solutions

[Next page](https://fenicsproject.discourse.group/t/unable-to-compute-the-curl-of-a-vector-field/14979.md?page=2)