I found that if I use mpirun I will get the write issue. The single core job will get rid of the issue. So there are something wrong with the mpi implementation, do you have any suggestion on how to fix the bug?
Meshio has no concept of MPI.
What you can do, is the following
from mpi4py import MPI
if MPI.COMM_WORLD.rank == 0:
# Do meshio code in here
MPI.COMM_WORLD.Barrier()
# Do dolfin stuff after this barrier
1 Like
I tried that. When I submitted the job, it seems the code did not really run. I even didnât see the first print of â-----------â as can be seen in the codes.
The reading in of the mesh should be outside of the loop, ie, the code after meshio.write(âŚ)
1 Like
That solved the issue. My last question is regarding the output. I use the command
File(âS1_V8_test1_results/displacement_test1.pvdâ) << u
to output the displacement. Then for parallel case, it seems each process will generate an .vtu file as can be seen in the screenshot. Is this normal or there is something I can improve for the output command? Thanks!
In general you are adviced to use other formats than pvd. For instance XDMFFile. In particular, if the function space is not a first order Lagrange space, you should use XDMFFile.write_checkpoint.
There are various topics on the subject here at discourse.
1 Like
Thank you so much! That solved the output issue. The specific code I used is: XDMFFile(âS1_V8_test1_results/displacement_test1.xdmfâ).write(u)
Hello,
How did you control the output format of the mesh? I got the same issue. In my case, I need to use the python interface of Gmsh to write the mesh. Do you have any idea how to specify the version II ASCII format using python script? Thank you.
Hello @dokken,
I try to convert the file test.msh (see below) into and xmdf file with the code you used above :
$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
7
2 29 "Left"
2 30 "Right"
2 31 "Surfaces"
3 25 "Iron"
3 26 "Silicon"
3 27 "Hoof"
3 28 "Bone"
$EndPhysicalNames
$ParametricNodes
12184
1 0 67.5 840 0 1
2 0 67.5 0 0 2
3 -7.654042494670956e-15 -67.5 0 0 3
...
12183 -15.56515802230132 -16.30129841836668 985 3 4
12184 -20.43860310699581 0.5179661288953525 985 3 4
$EndParametricNodes
$Elements
68727
1 15 2 0 1 1
2 15 2 0 2 2
...
68724 4 2 0 4 5855 15 5718 5598
68725 4 2 0 4 5855 5718 15 5578
68726 4 2 0 4 16 5718 15 5598
68727 4 2 0 4 16 15 5718 5578
$EndElements
However I get this error when trying to open my mesh file with meshio :
Without having the complete msh file, there is no way to give you any guidance. Please provide the msh file, or preferrably the geo file/python gmsh script that is used to generate the msh file.
Here are the step and geo files :
https://drive.google.com/drive/folders/1ImIVe21C7V2grRFhCmwj_mhUAyu7SnU2?usp=sharing
I did the procedure again, from the .step file et it seems to work now, I have a .xdmf.
Thanks for your reactivity !
Hello back dokken,
I have succeded in importing my .xdmf files :
I have succeded to apply boundary condition :
I would now like to create 4 subdomains corresponding to my 4 physical volumes labeled 25, 26, 27 and 28 in my .msh file. I have tried this :
But it does not works since fenics wants a mesh. If I use âmeshâ instead of âmaterialsâ in the last line, it doesnât recognise the different volumes. Could you please help me ?
Materials is already the mesh function you would use in your dx measure.
Please note that with only screenshots of snippets of code it is hard (and unlikely) that I or anyone else can resolve your issue.
Please put some effort into making a reproducible example or search the forum for similar topics, such as
Hello everyone, Iâm having trouble converting my .msh file to .xdmf. Iâd greatly appreciate it if anyone could offer some assistance. I created a rectangular shape in the GMSH interface and meshed it. Then, I exported it in version ASCÄ°Ä° 4 and saved it with the name ârectangler1.mshâ in my project folder in VS Code. However, when I run the code below, I get an error.
import meshio
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
# Algorithmic parameters
niternp = 20
niter = 80 # total number of iterations
pmax = 4
exponent_update_frequency = 4
tol_mass = 1e-4
thetamin = 0.001
# Problem parameters
thetamoy = 0.4
E = Constant(1)
nu = Constant(0.3)
lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/(2*(1+nu))
f = Constant((0, -1))
# Mesh
msh = meshio.read("rectangler1.msh")
tetra_cells = []
for cell in msh.cells:
if cell.type == "tetra":
if len(tetra_cells) == 0:
tetra_cells = cell.data
else:
tetra_cells = np.vstack([tetra_cells, cell.data])
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells})
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("mesh.xdmf", tetra_mesh)
mesh = Mesh("mesh.xdmf")
# Boundaries
def left(x, on_boundary):
return near(x[0], -2) and on_boundary
def load(x, on_boundary):
return near(x[0], 2) and near(x[1], 0.5, 0.05)
facets = MeshFunction("size_t", mesh, 1)
AutoSubDomain(load).mark(facets, 1)
ds = Measure("ds", subdomain_data=facets)
# Function space for density field
V0 = FunctionSpace(mesh, "DG", 0)
# Function space for displacement
V2 = VectorFunctionSpace(mesh, "CG", 2)
# Fixed boundary condtions
bc = DirichletBC(V2, Constant((0, 0)), left)
p = Constant(1) # SIMP penalty exponent
exponent_counter = 0 # exponent update counter
lagrange = Constant(1) # Lagrange multiplier for volume constraint
thetaold = Function(V0, name="Density")
thetaold.interpolate(Constant(thetamoy))
coeff = thetaold**p
theta = Function(V0)
volume = assemble(Constant(1.)*dx(domain=mesh))
avg_density_0 = assemble(thetaold*dx)/volume # initial average density
avg_density = 0.
def eps(v):
return sym(grad(v))
def sigma(v):
return coeff*(lamda*div(v)*Identity(2)+2*mu*eps(v))
def energy_density(u, v):
return inner(sigma(u), eps(v))
# Inhomogeneous elastic variational problem
u_ = TestFunction(V2)
du = TrialFunction(V2)
a = inner(sigma(u_), eps(du))*dx
L = dot(f, u_)*ds(1)
def local_project(v, V):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dx
b_proj = inner(v, v_)*dx
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
u = Function(V)
solver.solve_local_rhs(u)
return u
def update_theta():
theta.assign(local_project((p*coeff*energy_density(u, u)/lagrange)**(1/(p+1)), V0))
thetav = theta.vector().get_local()
theta.vector().set_local(np.maximum(np.minimum(1, thetav), thetamin))
theta.vector().apply("insert")
avg_density = assemble(theta*dx)/volume
return avg_density
def update_lagrange_multiplier(avg_density):
avg_density1 = avg_density
# Initial bracketing of Lagrange multiplier
if (avg_density1 < avg_density_0):
lagmin = float(lagrange)
while (avg_density < avg_density_0):
lagrange.assign(Constant(lagrange/2))
avg_density = update_theta()
lagmax = float(lagrange)
elif (avg_density1 > avg_density_0):
lagmax = float(lagrange)
while (avg_density > avg_density_0):
lagrange.assign(Constant(lagrange*2))
avg_density = update_theta()
lagmin = float(lagrange)
else:
lagmin = float(lagrange)
lagmax = float(lagrange)
# Dichotomy on Lagrange multiplier
inddico=0
while ((abs(1.-avg_density/avg_density_0)) > tol_mass):
lagrange.assign(Constant((lagmax+lagmin)/2))
avg_density = update_theta()
inddico += 1;
if (avg_density < avg_density_0):
lagmin = float(lagrange)
else:
lagmax = float(lagrange)
print(" Dichotomy iterations:", inddico)
def update_exponent(exponent_counter):
exponent_counter += 1
if (i < niternp):
p.assign(Constant(1))
elif (i >= niternp):
if i == niternp:
print("\n Starting penalized iterations\n")
if ((abs(compliance-old_compliance) < 0.01*compliance_history[0]) and
(exponent_counter > exponent_update_frequency) ):
# average gray level
gray_level = assemble((theta-thetamin)*(1.-theta)*dx)*4/volume
p.assign(Constant(min(float(p)*(1+0.3**(1.+gray_level/2)), pmax)))
exponent_counter = 0
print(" Updated SIMP exponent p = ", float(p))
return exponent_counter
u = Function(V2, name="Displacement")
old_compliance = 1e30
ffile = XDMFFile("topology_optimization.xdmf")
ffile.parameters["flush_output"]=True
ffile.parameters["functions_share_mesh"]=True
compliance_history = []
for i in range(niter):
solve(a == L, u, bc, solver_parameters={"linear_solver": "cg", "preconditioner": "hypre_amg"})
ffile.write(thetaold, i)
ffile.write(u, i)
compliance = assemble(action(L, u))
compliance_history.append(compliance)
print("Iteration {}: compliance =".format(i), compliance)
avg_density = update_theta()
update_lagrange_multiplier(avg_density)
exponent_counter = update_exponent(exponent_counter)
# Update theta field and compliance
thetaold.assign(theta)
old_compliance = compliance
plot(theta, cmap="bone_r")
plt.title("Final density")
plt.show();
plt.figure()
plt.plot(np.arange(1, niter+1), compliance_history)
ax = plt.gca()
ymax = ax.get_ylim()[1]
plt.plot([niternp, niternp], [0, ymax], "--k")
plt.annotate(r"$\leftarrow$ Penalized iterations $\rightarrow$", xy=[niternp+1, ymax*0.02], fontsize=14)
plt.xlabel("Number of iterations")
plt.ylabel("Compliance")
plt.title("Convergence history", fontsize=16)
plt.show();
zeynep@LAPTOP-3KN21JEE:~$ cd /home/zeynep ; /usr/bin/env /bin/python3 /home/zeynep/.vscode-server/extensions/ms-python.debugpy-2024.14.0-linux-x64/bundled/libs/debugpy/adapter/../../debugpy/launcher 42837 -- /home/zeynep/FEM/fem1.py
Traceback (most recent call last):
File "/home/zeynep/FEM/fem1.py", line 28, in <module>
msh = meshio.read("rectangler1.msh")
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3/dist-packages/meshio/_helpers.py", line 71, in read
return _read_file(Path(filename), file_format)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3/dist-packages/meshio/_helpers.py", line 90, in _read_file
raise ReadError(f"File {path} not found.")
meshio._exceptions.ReadError: File rectangler1.msh not found.
This isnât a DOLFIN issue, as the error message comes from MESHIO, and the fact that the file rectangler1.msh is not found, as clearly stated by the error message.
Please ensure that the file is in the same folder as you are calling your script from.
Secondly, as you havenât actually provided the mesh, the error would not be possible to reproduce by anyone else.
Finally, as the error happens early on in your code, you can massively simplify the code to:
import numpy as np
import meshio
msh = meshio.read("rectangler1.msh")
tetra_cells = []
cells) == 0:
tetra_cells = cell.data
else:
tetra_cells = np.vstack([tetra_cells, cell.data])
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells})
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("mesh.xdmf", tetra_mesh)
and the error would appear in line 3
Sorry, here is my mesh file.
Additionally, when I run the code just for the part you mentioned, I still encounter the same error.
Seemed like something when wrong when copy-pasting this. It should have been:
# Mesh
import meshio
import numpy as np
msh = meshio.read(filename="rectangler1.msh")
tetra_cells = []
for cell in msh.cells:
if cell.type == "tetra":
if len(tetra_cells) == 0:
tetra_cells = cell.data
else:
tetra_cells = np.vstack([tetra_cells, cell.data])
Running this, I get:
File "/root/shared/tmp-mesh/script.py", line 4, in <module>
msh = meshio.read(filename="rectangler1.msh")
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/dolfinx-env/lib/python3.12/site-packages/meshio/_helpers.py", line 71, in read
return _read_file(Path(filename), file_format)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/dolfinx-env/lib/python3.12/site-packages/meshio/_helpers.py", line 103, in _read_file
return reader_map[file_format](str(path))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/dolfinx-env/lib/python3.12/site-packages/meshio/gmsh/main.py", line 19, in read
mesh = read_buffer(f)
^^^^^^^^^^^^^^
File "/dolfinx-env/lib/python3.12/site-packages/meshio/gmsh/main.py", line 48, in read_buffer
return reader.read_buffer(f, is_ascii, data_size)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/dolfinx-env/lib/python3.12/site-packages/meshio/gmsh/_gmsh22.py", line 60, in read_buffer
has_additional_tag_data, cell_tags = _read_cells(
^^^^^^^^^^^^
File "/dolfinx-env/lib/python3.12/site-packages/meshio/gmsh/_gmsh22.py", line 136, in _read_cells
point_tags = np.asarray(point_tags, dtype=np.int32) - 1
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: int() argument must be a string, a bytes-like object or a real number, not 'NoneType'
which indicates that something is wrong with your .msh file.
As stated earlier, this is a GMSH or meshio question, and is better asked at:
Issues ¡ gmsh / gmsh ¡ GitLab (GMSH) or
Numerical Software (Meshio Discord)
1 Like