hi all,
im trying to create a mesh with my gmsh application installed on my pc and import it to xdmf
i was wondering if im creating a 2d model and extrude it to 3d, should i user tetra and triangle element or line and triangle ? (when i create the mesh in gmsh i do chose the option of 3d mesh)
im attaching my geo file and the mesh Dropbox
basically what i indented to do (as maybe some of you saw in my previous posts) is to create a rod with a cylinder instead of it (2 different materials).
im using an import that was posted before by Francesco Creating subdomains and boundary conditions on imported gmsh mesh
import os
import fenics as fe
from dolfin import *
import matplotlib.pyplot as plt
from ffc.fiatinterface import create_quadrature as cquad
import numpy as np
import meshio
import dolfin
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
#defining the young modulus
E1, E2, nu = 3e6, 3e6, 0.49
# Create mesh and define function space
mesh_name = "might_be_final"
# Then, read it in through meshio
meshio_mesh = meshio.read(f"{mesh_name}.msh")
# Extract mesh corresponding to each cell type, and write it out with meshio
def extract_mesh_by_cell_type(meshio_mesh, cell_type):
cells = meshio_mesh.get_cells_type(cell_type)
if len(cells) == 0:
return None
cell_data = meshio_mesh.get_cell_data("gmsh:physical", cell_type)
return meshio.Mesh(points=meshio_mesh.points, cells={cell_type: cells}, cell_data={"markers": [cell_data]})
tetra_mesh = extract_mesh_by_cell_type(meshio_mesh, "tetra")
if tetra_mesh:
meshio.write(f"{mesh_name}_tetra.xdmf", tetra_mesh)
triangle_mesh = extract_mesh_by_cell_type(meshio_mesh, "triangle")
if triangle_mesh:
meshio.write(f"{mesh_name}_triangle.xdmf", triangle_mesh)
del meshio_mesh
# Read back in the meshio mesh from XDMF in order to convert to a dolfin mesh
dolfin_mesh = dolfin.Mesh()
if triangle_mesh and os.path.exists(f"{mesh_name}_triagnle.xdmf"):
dolfin.XDMFFile(f"{mesh_name}_triangle.xdmf").read(dolfin_mesh)
elif tetra_mesh and os.path.exists(f"{mesh_name}_tetra.xdmf"):
dolfin.XDMFFile(f"{mesh_name}_tetra.xdmf").read(dolfin_mesh)
else:
raise FileNotFoundError("No suitable mesh file found to read.")
# Read back in the markers in order to attach them to the dolfin mesh
def read_back_in_markers(meshio_xdmf_file, dolfin_mesh):
dolfin_markers = dolfin.MeshValueCollection("size_t", dolfin_mesh)
dolfin.XDMFFile(meshio_xdmf_file).read(dolfin_markers, "markers")
dolfin_mesh_function = dolfin.cpp.mesh.MeshFunctionSizet(dolfin_mesh, dolfin_markers)
dolfin_mesh_function_array = dolfin_mesh_function.array()
defined_markers = np.array(list(set(dolfin_markers.values().values())), dtype=dolfin_mesh_function_array.dtype)
defined_entries = np.isin(dolfin_mesh_function_array, defined_markers)
dolfin_mesh_function_array[~defined_entries] = 0
dolfin_mesh_function.set_values(dolfin_mesh_function_array)
return dolfin_mesh_function
dolfin_subdomains = read_back_in_markers(f"{mesh_name}_tetra.xdmf", dolfin_mesh)
dolfin_boundaries = read_back_in_markers(f"{mesh_name}_triangle.xdmf", dolfin_mesh)
# Now export them to file
dolfin.XDMFFile(f"{mesh_name}.xdmf").write(dolfin_mesh)
dolfin.XDMFFile(f"{mesh_name}_subdomains.xdmf").write(dolfin_subdomains)
dolfin.XDMFFile(f"{mesh_name}_boundaries.xdmf").write(dolfin_boundaries)
# Please now open these three files in paraview and check that the labels are as expected
# You may remove *_tetra.* and *_triangle.*, they are not needed anymore
os.remove(f"{mesh_name}_tetra.xdmf")
os.remove(f"{mesh_name}_tetra.h5")
if triangle_mesh:
os.remove(f"{mesh_name}_triangle.xdmf")
os.remove(f"{mesh_name}_triangle.h5")
mesh_name = "might_be_final"
mesh = dolfin.Mesh()
dolfin.XDMFFile(f"{mesh_name}.xdmf").read(mesh)
materials = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim())
dolfin.XDMFFile(f"{mesh_name}_subdomains.xdmf").read(materials)
print(f"Subdomains: {np.unique(materials.array())}")
boundary_parts = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
dolfin.XDMFFile(f"{mesh_name}_boundaries.xdmf").read(boundary_parts)
print(f"Boundaries: {np.unique(boundary_parts.array())}")
tol = 1E-14
V = VectorFunctionSpace(mesh, "Lagrange", 1)
bc1 = DirichletBC(V.sub(0), Constant(0), boundary_parts, 19)
bc3 = DirichletBC(V.sub(1), Constant(0), boundary_parts, 21)
bc5 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 23)
bc6 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 24)
tDirBC = Expression(('1.0*time_'), time_=0.0, degree=0)
tDirBC4 = Expression(('((1/(time_+1))-1)*0.2'), time_=0.0, degree=0)
bc2 = DirichletBC(V.sub(0), tDirBC, boundary_parts, 20)
bc4 = DirichletBC(V.sub(1), tDirBC4, boundary_parts, 22)
bcs = [bc1, bc2, bc3, bc4, bc5, bc6]
subdomain_ids = [20, 22]
ds = Measure("ds", subdomain_data=boundary_parts, subdomain_id=(20, 22))
ds = ds(degree=4)
more then appreciate any advise/help
