Hello everyone,
I am using FEniCSX (fenics-dolfinx 0.9.0) with Linux (Ubuntu 24.04.2 LTS) which i installed with apt.
I create a custom made mesh using the code i found there Mesh generation
Now i modify the code from Mixed formulation for the Poisson equation.
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
div, exp, inner)
import matplotlib.pyplot as plt
import pyvista
import basix.ufl
import dolfinx
import ufl
#domain = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)
nodes = np.array([[1.0, 0.0], [2.0, 0.0], [3.0, 2.0], [1, 3]], dtype=np.float64)
connectivity = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int64)
c_el = ufl.Mesh(basix.ufl.element("CG", "triangle", 1, shape=(nodes.shape[1],)))
domain = dolfinx.mesh.create_mesh(MPI.COMM_SELF, connectivity, nodes, c_el)
k = 4
Q_el = element("BDM", domain.basix_cell(), k)
P_el = element("CG", domain.basix_cell(), k - 1)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(domain, V_el)
(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)
x = SpatialCoordinate(domain)
f = 10.0 * exp(-((x[0] - 0.5) * (x[0] - 0.5) + (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)
dx = Measure("dx", domain)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx
fdim = domain.topology.dim - 1
facets_top = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 1.0))
Q, _ = V.sub(0).collapse()
dofs_top = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_top)
def f1(x):
values = np.zeros((2, x.shape[1]))
values[1, :] = np.sin(5 * x[0])
return values
f_h1 = fem.Function(Q)
f_h1.interpolate(f1)
bc_top = fem.dirichletbc(f_h1, dofs_top, V.sub(0))
facets_bottom = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 0.0))
dofs_bottom = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_bottom)
def f2(x):
values = np.zeros((2, x.shape[1]))
values[1, :] = -np.sin(5 * x[0])
return values
f_h2 = fem.Function(Q)
f_h2.interpolate(f2)
bc_bottom = fem.dirichletbc(f_h2, dofs_bottom, V.sub(0))
bcs = [bc_top, bc_bottom]
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
"pc_factor_mat_solver_type": "mumps"})
try:
w_h = problem.solve()
except PETSc.Error as e: # type: ignore
if e.ierr == 92:
print("The required PETSc solver/preconditioner is not available. Exiting.")
print(e)
exit(0)
else:
raise e
sigma_h, u_h = w_h.split()
with io.XDMFFile(domain.comm, "out_mixed_poisson/u.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(u_h)
For k=1 is working but once i increase the value of k i get the following error
Traceback (most recent call last):
File "/home/vboxuser/Desktop/demo_mixed-poisson.py", line 82, in <module>
file.write_function(u_h)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/io/utils.py", line 252, in write_function
super().write_function(getattr(u, "_cpp_object", u), t, mesh_xpath)
RuntimeError: Degree of output Function must be same as mesh degree. Maybe the Function needs to be interpolated?
Even when i interpolate the solution into a Function space with the code
V_out = fem.functionspace(domain, element("CG", domain.basix_cell(), k))
u_h_interp = fem.Function(V_out)
u_h_interp.interpolate(u_h)
i get the same
Traceback (most recent call last):
File "/home/vboxuser/Desktop/demo_mixed-poisson.py", line 86, in <module>
file.write_function(u_h_interp)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/io/utils.py", line 252, in write_function
super().write_function(getattr(u, "_cpp_object", u), t, mesh_xpath)
RuntimeError: Degree of output Function must be same as mesh degree. Maybe the Function needs to be interpolated?
How can this be corrected?