hi,
i am new to dolfinx.
i have an external gmsh file which is in .msh format, i am trying to import it to dolfinx using meshio but i do not seem to see any mesh coming out. how do i know that it is imported and how can i see the mesh on jupyter lab.
please is there any guide or documentation on this.
Also, want to show flow only into and out of the mesh
i do not if something is wrong but i am trying the same codes that are on this [Defining subdomains for different materials — FEniCSx tutorial)
and getting this for the Subdomains defined from external mesh data
and i do not know what is missing and also for the first part the mesh created is not interactive it is static, how do i make it look interactive on jupyter lab
it is unable to open mesh.msh …am i supposed to download it from somewhere?
Try restarting your kernel, as the error you have there usually indicates that you have executed this cell before (that is why opencascade complains, as you have added this object once before). The mesh.msh gets generated when you call gmsh.write("mesh.msh").
1 Like
here is the link to the mesh https://shorturl.at/mxyJ6
so i am using the sames codes as the Test problem 1: Channel flow (Poiseuille flow) in the documentation but i got this error after running this code:
Define strain-rate tensor
def epsilon(u):
return sym(nabla_grad(u))
Define stress tensor
def sigma(u, p):
return 2muepsilon(u) - p*Identity(len(u))
Define the variational problem for the first step
p_n = Function(Q)
p_n.name = “p_n”
F1 = rho*dot((u - u_n) / k, v)dx
F1 += rhodot(dot(u_n, nabla_grad(u_n)), v)*dx
F1 += inner(sigma(U, p_n), epsilon(v))dx
F1 += dot(p_nn, v)ds - dot(munabla_grad(U)*n, v)*ds
F1 -= dot(f, v)*dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
the error:
Dimension mismatch in dot product.
ERROR:UFL:Dimension mismatch in dot product.
UFLException Traceback (most recent call last)
Cell In[32], line 16
14 F1 += inner(sigma(U, p_n), epsilon(v))dx
15 F1 += dot(p_nn, v)ds - dot(munabla_grad(U)*n, v)*ds
—> 16 F1 -= dot(f, v)*dx
17 a1 = form(lhs(F1))
18 L1 = form(rhs(F1))
File ~/mambaforge/envs/fenicsx-env/lib/python3.10/site-packages/ufl/operators.py:177, in dot(a, b)
175 if a.ufl_shape == () and b.ufl_shape == ():
176 return a * b
→ 177 return Dot(a, b)
File ~/mambaforge/envs/fenicsx-env/lib/python3.10/site-packages/ufl/tensoralgebra.py:193, in Dot.new(cls, a, b)
190 error("Dot product requires non-scalar arguments, "
191 “got arguments with ranks %d and %d.” % (ar, br))
192 if not (scalar or ash[-1] == bsh[0]):
→ 193 error(“Dimension mismatch in dot product.”)
195 # Simplification
196 if isinstance(a, Zero) or isinstance(b, Zero):
File ~/mambaforge/envs/fenicsx-env/lib/python3.10/site-packages/ufl/log.py:135, in Logger.error(self, *message)
133 “Write error message and raise an exception.”
134 self._log.error(*message)
→ 135 raise self._exception_type(self._format_raw(*message))
UFLException: Dimension mismatch in dot product.
which seems like the tutorial is for a 2D case and mine is 3D…is there a documentation that i can follow or what changes should i make to the code to be able to show this : with flow only into and out of the two holes, and no flow out the top or bottom (assuming effectively it’s a wall) of the mesh.
First of all, the link you provided does not work, and please do not use short-urls, as they are not transparent.
Secondly, please use 3x` encapsulation of code, i.e.
```python
import dolfinx
#add code here
```
Finally:
if your mesh is three-dimensional, f has to be 3D as well.
this is the link to the mesh https://drive.google.com/file/d/1PK5vIlw8lP4tE_A-9jmV9NYf7qxclKn1/view?usp=drive_link
and also which documentation explains f for 3D …i do not know how to do it
I assume you want to use Zero body forces?
If you look in the tutorial you are referring to:
https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code1.html
f is defined as
f = Constant(mesh, PETSc.ScalarType((0,0)))
Thus in 3D it should be
f = Constant(mesh, PETSc.ScalarType((0,0,0)))
thank you…i tried it and i got past that i am trying to achieve this flow. the diagram is attached below:
but with the codes i ran:
```python
import dolfinx
#from mpi4py import MPI
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh("layer.msh", MPI.COMM_WORLD, gdim=3)
from dolfinx.plot import create_vtk_mesh
import meshio
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]})
return out_mesh
proc = MPI.COMM_WORLD.rank
import numpy as np
import pyvista
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from ufl import (FacetNormal, FiniteElement, Identity,TestFunction, TrialFunction, VectorElement,
div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)
v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
def walls(x):
return np.logical_or(np.isclose(x[1],0), np.isclose(x[1],1))
wall_dofs = locate_dofs_geometrical(V, walls)
u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bc_noslip = dirichletbc(u_noslip, wall_dofs, V)
t = 0
T = 10
num_steps = 500
dt = T/num_steps
u_n = Function(V)
u_n.name = "u_n"
U = 0.5 * (u_n + u)
n = FacetNormal(mesh)
f = Constant(mesh, PETSc.ScalarType((0,0,0)))
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(1))
rho = Constant(mesh, PETSc.ScalarType(1))
def inflow(x):
return np.isclose(x[0], 0)
inflow_dofs = locate_dofs_geometrical(Q, inflow)
bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q)
def outflow(x):
return np.isclose(x[0], 1)
outflow_dofs = locate_dofs_geometrical(Q, outflow)
bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q)
bcu = [bc_noslip]
bcp = [bc_inflow, bc_outflow]
# Define strain-rate tensor
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))
# Define the variational problem for the first step
p_n = Function(Q)
p_n.name = "p_n"
F1 = rho*dot((u - u_n) / k, v)*dx
F1 += rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx
F1 += inner(sigma(U, p_n), epsilon(v))*dx
F1 += dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds
F1 -= dot(f, v)*dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = assemble_matrix(a1, bcs=bcu)
A1.assemble()
b1 = create_vector(L1)
# Define variational problem for step 2
u_ = Function(V)
a2 = form(dot(nabla_grad(p), nabla_grad(q))*dx)
L2 = form(dot(nabla_grad(p_n), nabla_grad(q))*dx - (rho/k)*div(u_)*q*dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)
# Define variational problem for step 3
p_ = Function(Q)
a3 = form(rho*dot(u, v)*dx)
L3 = form(rho*dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)
# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")
# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.BCGS)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")
# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)
pyvista.start_xvfb()
topology, cell_types, geometry = create_vtk_mesh(V)
values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
values[:, :len(u_n)] = u_n.x.array.real.reshape((geometry.shape[0], len(u_n)))
# Create a point cloud of glyphs
function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
function_grid["u"] = values
glyphs = function_grid.glyph(orient="u", factor=0.2)
# Create a pyvista-grid for the mesh
grid = pyvista.UnstructuredGrid(*create_vtk_mesh(mesh, mesh.topology.dim))
# Create plotter
plotter = pyvista.Plotter()
plotter.add_mesh(grid, style="wireframe", color="k")
plotter.add_mesh(glyphs)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
fig_as_array = plotter.screenshot("glyphs.png")
and i got this instead. i do not know if i am doing something wrong…is there something i needed to correct…i am trying to achieve the first picture i uploaded.
You haven’t meshed the interior of your domain, so how so how do you expect there to be any fluid flow in the interior?
okay sir…how do i go about doing that? is there any part of the code i am missing…
i am going through the tutorials. i do not know which one explains domain interior meshing…please can you guide me.
You need to mesh the entire domain you want to do simulations on, you cannot just mesh the boundaries. See for instance Using the GMSH Python API to generate complex meshes | Jørgen S. Dokken
i went through it but i really do not understand what and how to go about meshing the domain.
from this heading: How to create a 3D mesh with physical tags with the GMSH python API.
i believe i had created the 3D mesh right? so should i continue by running all the codes highlighted on this Using the GMSH Python API to generate complex meshes | Jørgen S. Dokken
From your figure it looks like you want to simulate fluid flow in the interior of your mesh:
I.e. in the area where you haven’t added any cells.
You would have to mesh the interior box if you want to simulate flow through it.
ok sir, so which documentation explains how to simulate the flow. i used the naiver stokes tutorial for poiseulle flow but it came out the earlier picture i uploaded and you suggested that i mesh the interior domain which you proposed i looked at this Using the GMSH Python API to generate complex meshes …
i am trying to understand which codes should i run? can you please guide me through it…
is i something that i am going to use gmsh to do or is there a particular codes on fenicsx to use to do that…i need help on how to go about it
You need to add the volume of the fluid to your mesh when creating it in GMSH. I do not have the time to guide you through this process, it should be rather straightforward as you have generated what seems like the enclosing surfaces of your volume.