2D Poisson with Vector Basis Functions

General
1

Hello,

I am new to Fenicsx and FEM, I studied on tutorials and previous discussions but I could not solve this problem, it is a simple question for you but I would be very happy if you answer it, it is important for me. I want to solve the 2D Poisson problem in the relevant mesh I created and I want to keep the stiffness matrix, right-hand side and solution as a matrix form. I wrote a code like the one below for this, but I still get the error in the rest of the code. Also, will my code enable 2D Poisson solving using Nedelec basis functions on the relevant mesh? Also, how can I add the 0 or another Dirichlet boundary condition to this scenario ? Thanks

import gmsh
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import mesh, fem, plot, io
from dolfinx.io import gmshio
import pyvista
from petsc4py import PETSc
import scipy.sparse

Initialize visualization

pyvista.start_xvfb()

def create_mesh_with_num_triangles(target_num_triangles):
if gmsh.isInitialized():
gmsh.finalize()
gmsh.initialize()
gmsh.model.add(“Mesh with target triangles”)
rectangle_tag = gmsh.model.occ.addRectangle(0, 0, 0, 1, 1)
gmsh.model.occ.synchronize()
average_edge_length = np.sqrt(2 / target_num_triangles)
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), average_edge_length)
gmsh.model.mesh.generate(2)
gmsh.model.addPhysicalGroup(2, [rectangle_tag], 1)
gmsh.write(“target_mesh.msh”)
gmsh.finalize()
MPI.COMM_WORLD.Barrier()
return gmshio.read_from_msh(“target_mesh.msh”, MPI.COMM_WORLD, gdim=2)[0]

def visualize_mesh(mesh):
topology, cell_types, geometry = plot.vtk_mesh(mesh, mesh.geometry.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
plotter.screenshot(“target_mesh.png”)

Create and visualize mesh

mesh = create_mesh_with_num_triangles(100)
visualize_mesh(mesh)

Setup function space and problem

V = fem.FunctionSpace(mesh, ufl.FiniteElement(“Nedelec 1st kind H(curl)”, mesh.ufl_cell(), 1))
sigma = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = fem.form(ufl.inner(ufl.curl(sigma), ufl.curl(v)) * ufl.dx)
L = fem.form(ufl.inner(ufl.as_vector((1, 1)), v) * ufl.dx)

Assembly system

A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a))
A.assemble()
b = fem.assemble_vector(L)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # Ensure correct assembly

Solver setup

solver = PETSc.KSP().create(MPI.COMM_WORLD)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
solver.setOperators(A)
u = fem.Function(V)

#Solve
solver.solve(b, u.vector)

Save solutions

with io.XDMFFile(MPI.COMM_WORLD, “solution.xdmf”, “w”) as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(u)

Error: ---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[16], line 56
54 A.assemble()
55 b = fem.assemble_vector(L)
—> 56 b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # Ensure correct assembly
58 # Solver setup
59 solver = PETSc.KSP().create(MPI.COMM_WORLD)

AttributeError: ‘Vector’ object has no attribute ‘ghostUpdate’

I used 0.7.2. Dolfinx version.

Thank you,

Chen

2

Please format your code properly, otherwise it may be very hard to follow. Use “```” to get the right formatting.

For the specific error you get, you’d probably want to assemble b as b = fem.petsc.assemble_vector(L) (note the additional petsc), otherwise the output vector b is not in PETSc format, and hence you cannot use the PETSc-specific method ghostUpdate.

3

Thank you sir, I wrote this code to solve curl (curlE) = J and Homogeneous Dirichlet (zero values on the boundaries) boudary conditions. Code is attached to the below. My question is, error of this solution is closer to 1, what is my error in the code. I am new in Fenicsx, if you help me I will be very happy. I searched it but looks correct. Normalized L2 norm of the error: 0.9530089726246179
The code:

"def create_mesh_with_num_triangles(target_num_triangles):
if gmsh.isInitialized():
gmsh.finalize()
gmsh.initialize()
gmsh.model.add(“Mesh with target triangles”)
rectangle_tag = gmsh.model.occ.addRectangle(0, 0, 0, 1, 1)
gmsh.model.occ.synchronize()
average_edge_length = np.sqrt(2 / target_num_triangles)
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), average_edge_length)
gmsh.model.mesh.generate(2)
gmsh.model.addPhysicalGroup(2, [rectangle_tag], 1)
gmsh.write(“target_mesh.msh”)
gmsh.finalize()
MPI.COMM_WORLD.Barrier()
return gmshio.read_from_msh(“target_mesh.msh”, MPI.COMM_WORLD, gdim=2)[0]

def visualize_mesh(mesh):
topology, cell_types, geometry = plot.vtk_mesh(mesh, mesh.geometry.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
plotter.screenshot(“target_mesh.png”)

Create and visualize mesh

target_num_triangles = 100
mesh = create_mesh_with_num_triangles(target_num_triangles)
visualize_mesh(mesh)

Define function space using Nedelec basis functions for H(curl) space

element = ufl.FiniteElement(“N1curl”, mesh.ufl_cell(), 1)
V = fem.FunctionSpace(mesh, element)

Define trial and test functions

E = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

Define source term (current density)

J = fem.Constant(mesh, PETSc.ScalarType((1.0, 0.0))) # Non-zero source term

Define variational problem

a = ufl.inner(ufl.curl(E), ufl.curl(v)) * ufl.dx
L = ufl.inner(J, v) * ufl.dx

Apply boundary conditions (homogeneous Dirichlet)

boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, lambda x: np.full(x.shape[1], True))
boundary_dofs = fem.locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets)

zero_vec = fem.Function(V)
with zero_vec.vector.localForm() as loc:
loc.set(0.0)
bc = fem.dirichletbc(zero_vec, boundary_dofs)

Assemble system

a_form = fem.form(a)
L_form = fem.form(L)
A = fem.petsc.assemble_matrix(a_form, bcs=[bc])
A.assemble()
b = fem.petsc.assemble_vector(L_form)
fem.petsc.apply_lifting(b, [a_form], [[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(b, [bc])

Solve the system

solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.CG)
solver.getPC().setType(PETSc.PC.Type.HYPRE)
solver.setTolerances(rtol=1e-8)

E_h = fem.Function(V)
solver.solve(b, E_h.vector)
E_h.x.scatter_forward()

Save the FEM solution vector

np.save(“solution_vector.npy”, E_h.vector.array)

Define the exact solution (since the exact solution for this problem is not trivial, we’ll use a known one)

x = ufl.SpatialCoordinate(mesh)
E_exact_expr = ufl.as_vector((ufl.sin(ufl.pi * x[0]), ufl.sin(ufl.pi * x[1])))

Interpolate the exact solution into the function space

E_exact = fem.Function(V)
E_exact_expr = fem.Expression(E_exact_expr, V.element.interpolation_points(), mesh.comm)
E_exact.interpolate(E_exact_expr)

Save the exact solution vector

np.save(“exact_solution_vector.npy”, E_exact.vector.array)

Compute the error in the original Nedelec space

error = fem.Function(V)
error.vector.array = E_h.vector.array - E_exact.vector.array

Compute the L2 norm of the error in the Nedelec space

L2_error_form = ufl.inner(error, error) * ufl.dx
error_L2_local = fem.assemble_scalar(fem.form(L2_error_form))
error_L2 = np.sqrt(mesh.comm.allreduce(error_L2_local, op=MPI.SUM))

Compute the L2 norm of the exact solution in the Nedelec space

L2_exact_form = ufl.inner(E_exact, E_exact) * ufl.dx
exact_L2_local = fem.assemble_scalar(fem.form(L2_exact_form))
exact_L2 = np.sqrt(mesh.comm.allreduce(exact_L2_local, op=MPI.SUM))

Normalize the L2 error

normalized_L2_error = error_L2 / exact_L2
print(f"Normalized L2 norm of the error in Nedelec space: {normalized_L2_error}")

Compute the maximum error at any degree of freedom in the Nedelec space

error_max = np.max(np.abs(E_exact.vector.array - E_h.vector.array))
print(f"Maximum error in Nedelec space: {error_max}“)”

4

The code is still not properly formatted and not runnable (e.g., imports are missing). No one will be able to help if you keep posting your code this way.

5

Sorry, for that, code is:

import gmsh
import numpy as np
import ufl
from mpi4py import MPI
import dolfinx
from dolfinx import mesh, fem, plot, io
from dolfinx.io import gmshio
import pyvista
from petsc4py import PETSc
import scipy.sparse
import gmsh
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import mesh, fem, plot, io
from petsc4py import PETSc
from dolfinx.fem.petsc import LinearProblem
from ufl import ds, dx, grad, inner
pyvista.start_xvfb()

def create_mesh_with_num_triangles(target_num_triangles):
    if gmsh.isInitialized():
        gmsh.finalize()
    gmsh.initialize()
    gmsh.model.add("Mesh with target triangles")
    rectangle_tag = gmsh.model.occ.addRectangle(0, 0, 0, 1, 1)
    gmsh.model.occ.synchronize()
    average_edge_length = np.sqrt(2 / target_num_triangles)
    gmsh.model.mesh.setSize(gmsh.model.getEntities(0), average_edge_length)
    gmsh.model.mesh.generate(2)
    gmsh.model.addPhysicalGroup(2, [rectangle_tag], 1)
    gmsh.write("target_mesh.msh")
    gmsh.finalize()
    MPI.COMM_WORLD.Barrier()
    return gmshio.read_from_msh("target_mesh.msh", MPI.COMM_WORLD, gdim=2)[0]

def visualize_mesh(mesh):
    topology, cell_types, geometry = plot.vtk_mesh(mesh, mesh.geometry.dim)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
    plotter = pyvista.Plotter()
    plotter.add_mesh(grid, show_edges=True)
    plotter.view_xy()
    if not pyvista.OFF_SCREEN:
        plotter.show()
    else:
        plotter.screenshot("target_mesh.png")

# Create and visualize mesh
target_num_triangles = 100
mesh = create_mesh_with_num_triangles(target_num_triangles)
visualize_mesh(mesh)

# Define function space using Nedelec basis functions for H(curl) space
element = ufl.FiniteElement("N1curl", mesh.ufl_cell(), 1)
V = fem.FunctionSpace(mesh, element)

# Define trial and test functions
E = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Define source term (current density)
J = fem.Constant(mesh, PETSc.ScalarType((1.0, 0.0)))  # Non-zero source term

# Define variational problem
a = ufl.inner(ufl.curl(E), ufl.curl(v)) * ufl.dx
L = ufl.inner(J, v) * ufl.dx

# Apply boundary conditions (homogeneous Dirichlet)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, lambda x: np.full(x.shape[1], True))
boundary_dofs = fem.locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets)

zero_vec = fem.Function(V)
with zero_vec.vector.localForm() as loc:
    loc.set(0.0)
bc = fem.dirichletbc(zero_vec, boundary_dofs)

# Assemble system
a_form = fem.form(a)
L_form = fem.form(L)
A = fem.petsc.assemble_matrix(a_form, bcs=[bc])
A.assemble()
b = fem.petsc.assemble_vector(L_form)
fem.petsc.apply_lifting(b, [a_form], [[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(b, [bc])

# Solve the system
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.CG)
solver.getPC().setType(PETSc.PC.Type.HYPRE)
solver.setTolerances(rtol=1e-8)

E_h = fem.Function(V)
solver.solve(b, E_h.vector)
E_h.x.scatter_forward()

# Save the FEM solution vector
np.save("solution_vector.npy", E_h.vector.array)

# Define the exact solution (since the exact solution for this problem is not trivial, we'll use a known one)
x = ufl.SpatialCoordinate(mesh)
E_exact_expr = ufl.as_vector((ufl.sin(ufl.pi * x[0]), ufl.sin(ufl.pi * x[1])))

# Interpolate the exact solution into the function space
E_exact = fem.Function(V)
E_exact_expr = fem.Expression(E_exact_expr, V.element.interpolation_points(), mesh.comm)
E_exact.interpolate(E_exact_expr)

# Save the exact solution vector
np.save("exact_solution_vector.npy", E_exact.vector.array)

# Compute the error in the original Nedelec space
error = fem.Function(V)
error.vector.array = E_h.vector.array - E_exact.vector.array

# Compute the L2 norm of the error in the Nedelec space
L2_error_form = ufl.inner(error, error) * ufl.dx
error_L2_local = fem.assemble_scalar(fem.form(L2_error_form))
error_L2 = np.sqrt(mesh.comm.allreduce(error_L2_local, op=MPI.SUM))

# Compute the L2 norm of the exact solution in the Nedelec space
L2_exact_form = ufl.inner(E_exact, E_exact) * ufl.dx
exact_L2_local = fem.assemble_scalar(fem.form(L2_exact_form))
exact_L2 = np.sqrt(mesh.comm.allreduce(exact_L2_local, op=MPI.SUM))

# Normalize the L2 error
normalized_L2_error = error_L2 / exact_L2
print(f"Normalized L2 norm of the error in Nedelec space: {normalized_L2_error}")

# Compute the maximum error at any degree of freedom in the Nedelec space
error_max = np.max(np.abs(E_exact.vector.array - E_h.vector.array))
print(f"Maximum error in Nedelec space: {error_max}")
6

I edited your latest post myself. Please go and have a look on how I did that, so that in future you don’t keep posting unformatted code.

General suggestions:

  • Are you really sure that J is the appropriate forcing term to have E_exact_expr as the exact solution? I haven’t done the math myself (I have never worked with curl/curl problems), but it’s worth double checking
  • If you are sure of what you wrote for J and E_exact_expr, what happens if you refine the mesh?