EM field from a meshed part

General
1

Dear all,

I am totally new to Fenics, so please excuse my low level questions …

I want to compute the electric field generated by of a meshed part which is at electrical potential of 100kV.
I have meshed th part with gmsh, and I would like to comute the electric field with dolfinX.

The meshed part is very simple (surface meshing):

image
image523×509 47.6 KB

So far, I have the following code:

import numpy as np
import matplotlib.pyplot as plt
import dolfinx
from mpi4py import MPI
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem import (FunctionSpace, dirichletbc, Function, Constant, locate_dofs_geometrical,
assemble_matrix, assemble_vector, set_bc)
from ufl import TrialFunction, TestFunction, inner, grad, dx
from petsc4py import PETSc

comm = MPI.COMM_WORLD

Load the mesh

with XDMFFile(comm, “generated_mesh_2D.xdmf”, “r”) as xdmf:
mesh = xdmf.read_mesh(name=“Grid”)

Define the function space

V = FunctionSpace(mesh, (“CG”, 1))

Define the boundary condition

u_D = Function(V)
u_D.interpolate(lambda x: np.full(x.shape[1], 100000)) # 100 kV

Boundary condition marker function

def boundary(x):
return np.full(x.shape[1], True)

dofs = locate_dofs_geometrical(V, boundary)
bc = dirichletbc(u_D, dofs)

Define the variational problem

u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
zero = Constant(mesh, 0.0) # Ensure the constant is properly initialized
L = zero * v * dx

Assemble and solve the system

A = assemble_matrix(a, bcs=[bc])
b = assemble_vector(L)
set_bc(b, [bc])

solver = PETSc.KSP().create(comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

phi = Function(V)
solver.solve(b, phi.vector)

Save the solution

with XDMFFile(MPI.COMM_WORLD, “electric_potential_solution.xdmf”, “w”) as file:
file.write_mesh(mesh)
file.write_function(phi)

But I get the following error:

AttributeError Traceback (most recent call last)
Cell In[38], line 40
37 L = zero * v * dx
39 # Assemble and solve the system
—> 40 A = assemble_matrix(a, bcs=[bc])
41 b = assemble_vector(L)
42 set_bc(b, [bc])

File ~/miniconda3/lib/python3.12/functools.py:909, in singledispatch..wrapper(*args, **kw)
905 if not args:
906 raise TypeError(f’{funcname} requires at least ’
907 ‘1 positional argument’)
→ 909 return dispatch(args[0].class)(*args, **kw)

File ~/miniconda3/lib/python3.12/site-packages/dolfinx/fem/assemble.py:243, in assemble_matrix(a, bcs, diagonal, constants, coeffs, block_mode)
219 “”“Assemble bilinear form into a matrix.
220
221 Args:
(…)
240
241 “””
242 bcs = if bcs is None else bcs
→ 243 A: la.MatrixCSR = create_matrix(a, block_mode)
244 _assemble_matrix_csr(A, a, bcs, diagonal, constants, coeffs)
245 return A

File ~/miniconda3/lib/python3.12/site-packages/dolfinx/fem/assemble.py:102, in create_matrix(a, block_mode)
92 def create_matrix(a: Form, block_mode: typing.Optional[la.BlockMode] = None) → la.MatrixCSR:
93 “”“Create a sparse matrix that is compatible with a given bilinear form.
94 Args:
95 a: Bilinear form to assemble.
(…)
100
101 “””
→ 102 sp = dolfinx.fem.create_sparsity_pattern(a)
103 sp.finalize()
104 if block_mode is not None:

File ~/miniconda3/lib/python3.12/site-packages/dolfinx/fem/init.py:37, in create_sparsity_pattern(a)
23 def create_sparsity_pattern(a: Form):
24 “”“Create a sparsity pattern from a bilinear form.
25
26 Args:
(…)
35
36 “””
—> 37 return _create_sparsity_pattern(a._cpp_object)

AttributeError: ‘Form’ object has no attribute ‘_cpp_object’

No idea of to deal with that.

First, I am not sure at all if this is the proper way to do it, and if by chance it is … how to solve this error message.

Could you please guide me during my first steps in dolfinx?

Many thanks in advance.

Best

2

Please format your code with 3x` syntax. Ie

```python
#add code here
```

and similarly for the error message.

Secondly,

Is not how you assemble in DOLFINx, i would suggest to have a look at
https://jsdokken.com/dolfinx-tutorial/chapter2/diffusion_code.html

3

Thanks for your reply, and sorry about the previous formating.

I have tried to adapt the code to deal properly with assemble, but unsuccessfully. Here is the new version of my code:

import numpy as np
import dolfinx
import dolfinx.fem as fem
import pyvista
from mpi4py import MPI
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem import (FunctionSpace, dirichletbc, Function, Constant, locate_dofs_geometrical,
                         assemble_matrix, create_vector, set_bc)
from ufl import TrialFunction, TestFunction, inner, grad, dx, dot
from petsc4py import PETSc

comm = MPI.COMM_WORLD

# Load the mesh
with XDMFFile(comm, "generated_mesh_2D.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")

# Define the scalar function space for potential
V = FunctionSpace(mesh, ("CG", 1))
# Define the vector function space for electric field
element = VectorElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, element)

# Define the boundary condition
u_D = Function(V)
u_D.interpolate(lambda x: np.full(x.shape[1], 100000))  # 100 kV

# Boundary condition marker function
def boundary(x):
    return np.full(x.shape[1], True)

# Locate degrees of freedom and apply boundary condition
boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim-1, boundary)
dofs = locate_dofs_geometrical(V, boundary)

# Apply Dirichlet boundary conditions
bc = dirichletbc(u_D, dofs)

# Define the variational problem for potential
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
zero = Constant(mesh, 0.0)
L = zero * v * dx

# Assemble the system
bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = create_vector(linear_form)
set_bc(b, [bc])

# Solve the system
solver = PETSc.KSP().create(comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

# Solve for potential
phi = Function(V)
solver.solve(b, phi.vector)

# Compute the electric field as the negative gradient of the potential
E = Function(W)
E.interpolate(lambda x: -grad(phi)(x))  # Electric field is negative gradient of potential

# Save the solution
with XDMFFile(comm, "electric_potential_solution.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(E)

# Convert to PyVista mesh for visualization
pvd_file = "electric_field.pvd"
with dolfinx.io.VTXWriter(comm, mesh, W, pvd_file, "w") as writer:
    writer.write(E)

# Read and plot using PyVista
grid = pyvista.read(pvd_file)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.add_arrows(grid.points, grid.point_data["Electric Field"], mag=1)
plotter.show()

I guess I still don’t do it properly.

Could you help me solving this issue?

Also, is it the proper way to do it for the plotting ?

Many thanks

4

You are still not creating the correct rhs.
Here you need to assemble the vector, then apply lifting, then set bc, as shown in:
http://jsdokken.com/FEniCS23-tutorial/src/lifting.html#lifting

For the remainder of your questions, you are not adding enough context for me to help you.

What do yo mean by this

5

Hello,

sorry I can’t make it run…`

Again, my idea here is the compute the electric field from a meshed imported part. This part is fixed to an electrical potential of 100kV → I just want to compute the 3D electric field (as a first step…).

I am still doing something wrong. The command “A_lifting.assemble()” does not work (AttributeError: ‘MatrixCSR’ object has no attribute ‘assemble’ ), so I guess I am doing something wrong before that.

Here is the updated code:

import numpy as np
import dolfinx
import dolfinx.fem as fem
import pyvista
from mpi4py import MPI
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem import (FunctionSpace, dirichletbc, Function, Constant, locate_dofs_geometrical,
                         assemble_matrix, assemble_vector, set_bc, apply_lifting)
from ufl import TrialFunction, TestFunction, inner, grad, dx, ds
from petsc4py import PETSc
import ufl

comm = MPI.COMM_WORLD

# Load the mesh
with XDMFFile(comm, "generated_mesh_2D.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")

# Define the scalar function space for potential
V = FunctionSpace(mesh, ("CG", 1))

# Define the boundary condition
u_D = Function(V)
u_D.interpolate(lambda x: np.full(x.shape[1], 100000))  # 100 kV

# Boundary condition marker function
def boundary(x):
    return np.full(x.shape[1], True)

# Locate degrees of freedom based on the boundary marker
boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim-1, boundary)
dofs = locate_dofs_geometrical(V, boundary)
# Apply Dirichlet boundary conditions
bc = dirichletbc(u_D, dofs)

# Define the variational problem for potential
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
f = Constant(mesh, 0.0)
L = f * v * dx

# Assemble the system
a_compiled  = fem.form(a)
L_compiled  = fem.form(L)

A_lifting  = assemble_matrix(a_compiled , bcs=[bc])
#A_lifting.assemble()
b_lifting  = assemble_vector(L_compiled )
apply_lifting(b_lifting, [a_compiled], bcs=[bc])
b_lifting.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b_lifting , [bc])

# Solve the system
solver = PETSc.KSP().create(comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

# Solve for potential
phi = Function(V)
solver.solve(b, phi.vector)

# Define the vector function space for electric field
element = fem.VectorElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, element)

# Compute the electric field as the negative gradient of the potential
E = Function(W)
E.interpolate(lambda x: -grad(phi)(x))  # Electric field is negative gradient of potential

# Save the solution
with XDMFFile(comm, "electric_potential_solution.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(phi)
    file.write_function(E)

# Convert to PyVista mesh for visualization
pvd_file = "electric_field.pvd"
with dolfinx.io.VTXWriter(comm, mesh, W, pvd_file, "w") as writer:
    writer.write(E)

# Read and plot using PyVista
grid = pyvista.read(pvd_file)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.add_arrows(grid.points, grid.point_vectors["E"], scale=1)
plotter.show()

could you please tell what I am doing wrong ? Unfortunately, I can not upload my mesh file here, but I guess any mesh file would work…

Thank you

6

Split this into


from dolfinx.fem import (FunctionSpace, dirichletbc, Function, Constant, locate_dofs_geometrical)
from dolfinx.fem.petsc import (assemble_matrix, assemble_vector, set_bc, apply_lifting)
7

Thanks for your reply.

Thanks to your help I made a small step forward… but I am now stuck with the E field computation.
I am pretty sure that I am doing it the wrong way.

The part where I am doing it is here:

element = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, element)
# Compute the electric field as the negative gradient of the potential
E = Function(W)
E_expr = -grad(phi)
# Save the solution
with XDMFFile(comm, "electric_potential_solution.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(phi)
    file.write_function(E_expr)

The code crashes as soon as I try to “file.write_function(E_expr)”, so I concluded I compute E the wrong way.

Any idea how to make it the proper way ?

Full code is here:

import numpy as np
import dolfinx
import dolfinx.fem as fem
import pyvista
from mpi4py import MPI
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem import (VectorFunctionSpace, FunctionSpace, dirichletbc, Function, Constant, locate_dofs_geometrical)
from dolfinx.fem.petsc import (assemble_matrix, assemble_vector, set_bc, apply_lifting)
from ufl import TrialFunction, TestFunction, inner, grad, dx, ds
from petsc4py import PETSc
import ufl
from ufl import VectorElement 
from ufl import FiniteElement 
from dolfinx.io import VTXWriter

comm = MPI.COMM_WORLD

# Load the mesh
with XDMFFile(comm, "generated_mesh_3D_Surface.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")

# Define the scalar function space for potential
V = FunctionSpace(mesh, ("CG", 1))

# Define the boundary condition
u_D = Function(V)
u_D.interpolate(lambda x: np.full(x.shape[1], 100000))  # 100 kV

# Boundary condition marker function
def boundary(x):
    return np.full(x.shape[1], True)

# Locate degrees of freedom based on the boundary marker
boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim-1, boundary)
dofs = locate_dofs_geometrical(V, boundary)
# Apply Dirichlet boundary conditions
bc = dirichletbc(u_D, dofs)

# Define the variational problem for potential
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
f = Constant(mesh, 0.0)
L = f * v * dx

# Assemble the system
a_compiled  = fem.form(a)
L_compiled  = fem.form(L)

#Lifting:
A_lifting  = assemble_matrix(a_compiled , bcs=[bc])
A_lifting.assemble()
b_lifting  = assemble_vector(L_compiled )
apply_lifting(b_lifting, [a_compiled], bcs=[[bc]])
b_lifting.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b_lifting , [bc])
print(f"Matrix is symmetric after lifting assembly: {A_lifting.isSymmetric(1e-5)}")

# Solve the system
solver = PETSc.KSP().create(comm)
solver.setOperators(A_lifting)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

# Solve for potential
phi = Function(V)
solver.solve(b_lifting, phi.vector)

# Define the vector function space for electric field
#element = VectorElement("CG", mesh.ufl_cell(), 1)
element = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, element)

# Compute the electric field as the negative gradient of the potential
E = Function(W)
#E.interpolate(lambda x: -grad(phi)(x))  # Electric field is negative gradient of potential
E_expr = -grad(phi)  # A more UFL-compliant expression
#E= fem.project(E_expr, W)  # Projecting the expression onto the function space W

print(f"E solved.")

# Save the solution
with XDMFFile(comm, "electric_potential_solution.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(phi)
    file.write_function(E_expr)

print(f"Xdmf file written.")

# Convert to PyVista mesh for visualization
pvd_file = "electric_field.pvd"
#with dolfinx.io.VTXWriter(comm, mesh, W, pvd_file, "w") as writer:
#    writer.write(E)
writer = VTXWriter(mesh, E_expr, pvd_file)  # File mode is optional and defaults usually to "w"
with writer:
    writer.write(E_expr)


print(f"PyVista convertion done.")

# Read and plot using PyVista
grid = pyvista.read(pvd_file)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.add_arrows(grid.points, grid.point_vectors["E"], scale=1)
plotter.show()

Thanks a lot for your help.

Best

8

Please see for instance: Elasticity using algebraic multigrid — DOLFINx 0.7.3 documentation