@dokken Maybe you can help me out. I have created simple square geometry in Gmsh and converted the .msh file into xdmf. When i import it in fenics using your code I’m getting different results compare to the simple command "UnitSquareMesh(30,30). I’m trying tonsolve the simple linear elastic problem. Maybe you can help me out. How i can import the external 2d mesh into fenics? is it the correct code to import 2d mesh? mesh = Mesh() with XDMFFile(“lavender.xdmf”) as infile: infile.read(mesh) mvc = MeshValueCollection(“size_t”, mesh, 1) with XDMFFile(“bound.xdmf”) as infile: infile.read(mvc, “SurfaceRegions”) mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

here you go to generate the square mesh: //+ Point(1) = {0, 0, 0, 1.0}; //+ Point(2) = {30, 0, 0, 1.0}; //+ Point(3) = {30, 30, 0, 1.0}; //+ Point(4) = {0, 30, 0, 1.0}; //+ Line(1) = {1, 2}; //+ Line(2) = {2, 3}; //+ Line(3) = {3, 4}; //+ Line(4) = {4, 1}; //+ Curve Loop(1) = {3, 4, 1, 2}; //+ Plane Surface(1) = {1}; //+ Physical Curve(5) = {1}; //+ Physical Curve(6) = {2}; //+ Physical Curve(7) = {3}; //+ Physical Curve(8) = {4}; //+ Physical Surface(9) = {1}; code to generate the mesh in xdmf import meshio mesh_from_file = meshio.read(“mesh.msh”) import numpy 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 line_mesh = create_mesh(mesh_from_file, “line”, prune_z=True) meshio.write(“bound.xdmf”, line_mesh) triangle_mesh = create_mesh(mesh_from_file, “triangle”, prune_z=True) meshio.write(“lavender.xdmf”, triangle_mesh) Fenics code: from dolfin import * from fenics import * from ufl import nabla_grad from ufl import nabla_div import numpy as np mesh=Mesh() with XDMFFile(‘lavender.xdmf’) as infile: infile.read(mesh) mvc = MeshValueCollection(“size_t”, mesh, 2) with XDMFFile(‘bound.xdmf’) as infile: infile.read(mvc,“name_to_read”) sub = cpp.mesh.MeshFunctionSizet(mesh, mvc) # Strain function def epsilon(u): return 0.5*(grad(u) + grad(u).T) #return sym(grad(u)) # Stress function def sigma(u): return lmbda_*div(u)*Identity(2) + 2*mu*epsilon(u) # Define material properties E = Constant(100000) nu = Constant(0.3) mu = E/(2*(1+nu)) lmbda_ = E*nu/((1+nu)*(1-2*nu)) flag_quad=2 P2 = VectorElement("Lagrange", mesh.ufl_cell(), flag_quad) TH = P2 W = FunctionSpace(mesh, TH) # Define traction on the boundary and body forces T = Constant((0.0, 0.0)) B = Constant((0.0, 0.0)) u = Function(W) du = TrialFunction(W) v = TestFunction(W) # Boundary Conditions def clamped_boundary(x, on_boundary): return on_boundary and x[1] < tol bc1 = DirichletBC(W, Constant((0, 0)), clamped_boundary) tol=1E-5 top = CompiledSubDomain('on_boundary && near(x[1], 1, tol)', tol=1E-5) applied_disp = 0.5 bc2 = DirichletBC(W.sub(1), Constant((applied_disp/2.0)), top) bcs = [bc1,bc2] a = inner(sigma(du), epsilon(v))*dx l = dot(B, v)*dx + inner(T, v)*ds(1) A_ass, L_ass = assemble_system(a, l, bcs) solve(A_ass, u.vector(), L_ass) plot(u, title="Displacement", mode="displacement") to_print = True file_results = XDMFFile("+-elasticity_reet.xdmf") file_results.parameters["flush_output"] = True file_results.parameters["functions_share_mesh"] = True file_results.write(u,0.)

Please edit your posts and make sure They use 3x encapsulation surrounding codes or files, i.e. ``` //+ Point(….) ``` and ```python from dolfin import * #…. ```

here you go to generate the square mesh: //+ Point(1) = {0, 0, 0, 1.0}; //+ Point(2) = {30, 0, 0, 1.0}; //+ Point(3) = {30, 30, 0, 1.0}; //+ Point(4) = {0, 30, 0, 1.0}; //+ Line(1) = {1, 2}; //+ Line(2) = {2, 3}; //+ Line(3) = {3, 4}; //+ Line(4) = {4, 1}; //+ Curve Loop(1) = {3, 4, 1, 2}; //+ Plane Surface(1) = {1}; //+ Physical Curve(5) = {1}; //+ Physical Curve(6) = {2}; //+ Physical Curve(7) = {3}; //+ Physical Curve(8) = {4}; //+ Physical Surface(9) = {1}; code to generate the mesh in xdmf import meshio mesh_from_file = meshio.read(“mesh.msh”) import numpy 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 line_mesh = create_mesh(mesh_from_file, “line”, prune_z=True) meshio.write(“bound.xdmf”, line_mesh) triangle_mesh = create_mesh(mesh_from_file, “triangle”, prune_z=True) meshio.write(“lavender.xdmf”, triangle_mesh) Fenics code: from dolfin import * from fenics import * from ufl import nabla_grad from ufl import nabla_div import numpy as np mesh=Mesh() with XDMFFile(‘lavender.xdmf’) as infile: infile.read(mesh) mvc = MeshValueCollection(“size_t”, mesh, 2) with XDMFFile(‘bound.xdmf’) as infile: infile.read(mvc,“name_to_read”) sub = cpp.mesh.MeshFunctionSizet(mesh, mvc) # Strain function def epsilon(u): return 0.5*(grad(u) + grad(u).T) #return sym(grad(u)) # Stress function def sigma(u): return lmbda_*div(u)*Identity(2) + 2*mu*epsilon(u) # Define material properties E = Constant(100000) nu = Constant(0.3) mu = E/(2*(1+nu)) lmbda_ = E*nu/((1+nu)*(1-2*nu)) flag_quad=2 P2 = VectorElement("Lagrange", mesh.ufl_cell(), flag_quad) TH = P2 W = FunctionSpace(mesh, TH) # Define traction on the boundary and body forces T = Constant((0.0, 0.0)) B = Constant((0.0, 0.0)) u = Function(W) du = TrialFunction(W) v = TestFunction(W) # Boundary Conditions def clamped_boundary(x, on_boundary): return on_boundary and x[1] < tol bc1 = DirichletBC(W, Constant((0, 0)), clamped_boundary) tol=1E-5 top = CompiledSubDomain('on_boundary && near(x[1], 1, tol)', tol=1E-5) applied_disp = 0.5 bc2 = DirichletBC(W.sub(1), Constant((applied_disp/2.0)), top) bcs = [bc1,bc2] a = inner(sigma(du), epsilon(v))*dx l = dot(B, v)*dx + inner(T, v)*ds(1) A_ass, L_ass = assemble_system(a, l, bcs) solve(A_ass, u.vector(), L_ass) plot(u, title="Displacement", mode="displacement") to_print = True file_results = XDMFFile("+-elasticity_reet.xdmf") file_results.parameters["flush_output"] = True file_results.parameters["functions_share_mesh"] = True file_results.write(u,0.)

[image] taytay33: top = CompiledSubDomain('on_boundary && near(x[1], 1, tol)', tol=1E-5) It is quite obvious as you have hardcoded the top boundary condition, which assumes that your mesh is a unit square, while what you actually have in your mesh is a [0,30]\times[0,30] mesh. Thus using top = CompiledSubDomain('on_boundary && near(x[1], 30, tol)', tol=1E-5) Here is a cleaned up version of your code: import matplotlib.pyplot as plt from dolfin import * import meshio mesh_from_file = meshio.read("mesh.msh") 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 line_mesh = create_mesh(mesh_from_file, "line", prune_z=True) meshio.write("bound.xdmf", line_mesh) triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True) meshio.write("lavender.xdmf", triangle_mesh) mesh = Mesh() with XDMFFile("lavender.xdmf") as infile: infile.read(mesh) mvc = MeshValueCollection("size_t", mesh, 2) with XDMFFile("bound.xdmf") as infile: infile.read(mvc, "name_to_read") sub = cpp.mesh.MeshFunctionSizet(mesh, mvc) # Strain function def epsilon(u): return 0.5*(grad(u) + grad(u).T) # return sym(grad(u)) # Stress function def sigma(u): return lmbda_*div(u)*Identity(2) + 2*mu*epsilon(u) # Define material properties E = Constant(100000) nu = Constant(0.3) mu = E/(2*(1+nu)) lmbda_ = E*nu/((1+nu)*(1-2*nu)) flag_quad = 2 P2 = VectorElement("Lagrange", mesh.ufl_cell(), flag_quad) TH = P2 W = FunctionSpace(mesh, TH) # Define traction on the boundary and body forces T = Constant((0.0, 0.0)) B = Constant((0.0, 0.0)) u = Function(W) du = TrialFunction(W) v = TestFunction(W) # Boundary Conditions def clamped_boundary(x, on_boundary): return on_boundary and x[1] < tol bc1 = DirichletBC(W, Constant((0, 0)), clamped_boundary) tol = 1E-5 top = CompiledSubDomain('on_boundary && near(x[1], 30, tol)', tol=1E-5) applied_disp = 0.5 bc2 = DirichletBC(W.sub(1), Constant((applied_disp/2.0)), top) bcs = [bc1, bc2] a = inner(sigma(du), epsilon(v))*dx l = dot(B, v)*dx + inner(T, v)*ds(1) A_ass, L_ass = assemble_system(a, l, bcs) solve(A_ass, u.vector(), L_ass) plot(u, title="Displacement", mode="displacement") plt.savefig("u.png") file_results = XDMFFile("+-elasticity_reet.xdmf") file_results.write(u, 0.) file_results.close() Note that I’ve made sure that the indentation of all functions are correct, which wasn’t the case in your code [image] taytay33: 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 Also note that you could use the marked boundaries from the gmsh mesh (“bound.xdmf”) to set boundary conditions.

@dokken thanks for your kind help. I really appreciate that. A little modification. I have this image domain in which i have 2 different materials the circles in domain have different material. how can I give different values of E to different domains in fenics? The applied boundary conditions are same as above: circles have E = 1e6 and othe part has E= 2e3 def clamped_boundary(x, on_boundary): return on_boundary and x[1] < tol bc1 = DirichletBC(W, Constant((0, 0)), clamped_boundary) tol = 1E-5 top = CompiledSubDomain('on_boundary && near(x[1], 30, tol)', tol=1E-5) applied_disp = 0.5 bc2 = DirichletBC(W.sub(1), Constant((applied_disp/2.0)), top) bcs = [bc1, bc2] Everything is same as above i have attached the geometry below: [different material1]

One question do i have to change this xdmf file? import matplotlib.pyplot as plt from dolfin import * import meshio mesh_from_file = meshio.read("mesh.msh") 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 line_mesh = create_mesh(mesh_from_file, "line", prune_z=True) meshio.write("bound.xdmf", line_mesh) triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True) meshio.write("lavender.xdmf", triangle_mesh)

[image] taytay33: 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]}) As long as you have Added physical tags in Gmsh, They are now saved in the xdmf files.

Let’s suppose the image given above is meshed using GMSH. The circles in images are extracted using OpenCv Contours detection. And I add physical tags to only boundary than whats the procedure to detect subdomains. In simple I use Gmsh for meshing the Image which has different domains how will fenics detect the subdoamins?

You can loop over the read in facets, find the connecting cells on each side and mark them. However, I would strongly advice you to marker the surfaces/subdomains (cell volumes) in Gmsh/your mesh generator, as it is alot more expensive to do recursive searches through neighbouring cells if you do not have specific geometrical information.

How to import xdmf mesh generated by gmsh into fenics

General

dokken December 29, 2022, 7:50am 12

See How to extract and use physical labels from GMSH surfaces on FEniCS (2D mesh) - #5 by dokken
and
How to define different materials / importet 3D geometry (gmsh) - #2 by dokken

Topic		Replies	Views
Assigning different physical groups to different color in gmsh and mark them as different material in fenics General	34	3635	June 6, 2023
Want to import already defined physical groups of Microstructure in Gmsh to FEnics mesh mesh	40	1766	January 13, 2024
Need help converting GMSH to FEniCS installation	61	7312	September 28, 2022
How to define different materials / importet 3D geometry (gmsh) Uncategorized	14	5512	November 4, 2020
Extracting mesh regions for gmsh mesh mesh	24	5554	December 5, 2023

How to import xdmf mesh generated by gmsh into fenics

Related topics