Hello everyone, I’m encountering an issue while attempting to plot my solutions. Both velocity and pressure appear to be null. I’m uncertain whether this issue stems from the geometry imported from Gmsh or other factors.
Below, you’ll find my code, the Gmsh file containing the geometry, and the velocity and pressure graphs visualized in ParaView
Thank you.
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import dolfinx.fem.petsc
import dolfinx.io
from basix.ufl import element, mixed_element
from dolfinx import mesh, fem
from dolfinx.fem import (Function, dirichletbc, form, functionspace,
locate_dofs_topological)
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner
import matplotlib.pyplot as plt
from dolfinx import plot
import pyvista
from dolfinx.plot import vtk_mesh
import warnings
warnings.filterwarnings("ignore")
import meshio
from dolfinx.io import XDMFFile
import numpy as np
import pyvista
from dolfinx.io import VTKFile
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data.astype(np.int32)]})
return out_mesh
msh = meshio.read("geo_cylindrique.msh")
for cell in msh.cells:
if cell.type == "triangle":
triangle_cells = cell.data
elif cell.type == "tetra":
tetra_cells = cell.data
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
elif key == "tetra":
tetra_data = msh.cell_data_dict["gmsh:physical"][key]
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells})
triangle_mesh =meshio.Mesh(points=msh.points,
cells=[("triangle", triangle_cells)],
cell_data={"name_to_read":[triangle_data]})
meshio.write("msh.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh)
with XDMFFile(MPI.COMM_WORLD, "msh.xdmf", "r") as xdmf:
msh = xdmf.read_mesh(name="Grid")
msh.topology.create_connectivity(msh.topology.dim, msh.topology.dim - 1)
with XDMFFile(MPI.COMM_WORLD, "mf.xdmf", "r") as xdmf:
ft = xdmf.read_meshtags(msh, name="Grid")
def u_expression(x):
R=0.5
Q=49.5
u_max=(3*Q)/(4*np.pi*R)
#(1-((x[1]**2+x[2]**2)/(R**2)))*(u_max)
return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1]),np.zeros(x.shape[1])))
def Dirichlet_inlet(msh,W,W0):
u_in=Function(W0)
u_in.interpolate(u_expression)
where=ft.indices[ft.values==5]
dofs = locate_dofs_topological((W.sub(0), W0), 2, where)
bc0 = dirichletbc(u_in, dofs, W.sub(0))
return bc0
def Dirichlet_walls(msh,W,W0):
u_wall= Function(W0)
where2=ft.indices[ft.values==7]
dofs1 = locate_dofs_topological((W.sub(0), W0), 2, where2)
bc1 = dirichletbc(u_wall, dofs1, W.sub(0))
return bc1
def mixed_direct(TH, msh):
W = functionspace(msh, TH)
W0, _ = W.sub(0).collapse()
W1, W1_to_W = W.sub(1).collapse()
bc0=Dirichlet_inlet(msh,W,W0)
bc1=Dirichlet_walls(msh,W,W0)
bcs = [bc0,bc1]
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Function(W0)
a = form((inner(grad(u), grad(v)) - inner(p, div(v)) + inner(div(u), q)) * dx)
L = form(inner(f, v) * dx)
A = fem.petsc.assemble_matrix(a, bcs=bcs)
A.assemble()
b = fem.petsc.assemble_vector(L)
fem.petsc.apply_lifting(b, [a], bcs=[bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(b, bcs)
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
U = Function(W)
try:
ksp.solve(b, U.vector)
except PETSc.Error as e:
if e.ierr == 92:
print("The required PETSc solver/preconditioner is not available. Exiting.")
print(e)
exit(0)
else:
raise e
ksp.solve(b, U.vector)
u, p = U.sub(0).collapse(), U.sub(1).collapse()
u.x.scatter_forward()
p.x.scatter_forward()
return u, p
if __name__ == "__main__":
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)
V = fem.functionspace(msh, P2)
TH = mixed_element([P2, P1])
u_h,p_h=mixed_direct(TH,msh)
with VTKFile(msh.comm, "3d_vitesse.pvd", "w") as file:
file.write_mesh(msh)
file.write_function([u_h._cpp_object])
with VTKFile(msh.comm, "3d_sol_pression.pvd", "w") as file:
file.write_mesh(msh)
file.write_function([p_h._cpp_object])
The gmsh script is belllow:
// Gmsh project created on Tue May 14 12:38:03 2024
SetFactory("OpenCASCADE");
//+
Cylinder(1) = {0, 0, 0, 4, 0, 0, 0.5, 2*Pi};
//+
Physical Volume(4) = {1};
//+
Physical Surface("inlet", 5) = {3};
//+
Physical Surface("outlet", 6) = {2};
//+
Physical Surface("walls", 7) = {1};
the solution pressure solution:
the velocity solution in P1 :
