PhsysicsSimGmsh.py

import gmsh
from mpi4py import MPI
from dolfinx.io import gmshio

# Initialize Gmsh and open the STL file
gmsh.initialize()
gmsh.model.add("Mesh")
gmsh.merge("DemoTO.stl")
#gmsh.model.mesh.import_stl()
gmsh.model.occ.synchronize()

# Generate 3D mesh
gmsh.model.mesh.generate(3)

# Export mesh as a .msh file
gmsh.write("mesh.msh")
gmsh.finalize()

# Read the mesh using dolfinx
from dolfinx.io import XDMFFile
from dolfinx.mesh import Mesh

with XDMFFile(MPI.COMM_WORLD, "mesh.msh", "r") as xdmf:
    mesh = xdmf.read_mesh()

from dolfinx import fem
from dolfinx.fem import Function, VectorFunctionSpace, Constant, dirichletbc
from dolfinx.mesh import locate_entities_boundary
from mpi4py import MPI
import ufl

# Define function space (Displacement)
V = VectorFunctionSpace(mesh, ("CG", 1))

# Define boundary condition
tol = 1E-14

# Apply zero displacement on a specific region (supports)
def boundary(x):
    return x[0] < tol  # Adjust based on your specific boundary conditions

# Locate boundary facets
boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, boundary)

# Create a zero displacement boundary condition
bc_value = fem.Function(V)
with bc_value.vector.localForm() as loc:
    loc.set(0.0)
bc = dirichletbc(bc_value, fem.locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets))

# Define force (load) applied to another region
force = fem.Constant(mesh, (0.0, 0.0, -10.0))  # Example force in the z-direction

# Define Lame's constants (material properties)
mu = fem.Constant(mesh, 10.0)  # Shear modulus
lambda_ = fem.Constant(mesh, 100.0)  # First Lame parameter

# Define strain and stress functions
def epsilon(u):
    return ufl.sym(ufl.grad(u))

def sigma(u):
    return lambda_ * ufl.div(u) * ufl.Identity(3) + 2 * mu * epsilon(u)

# Define variational problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(force, v) * ufl.ds

# Solve the problem
u_sol = fem.Function(V)
problem = fem.petsc.LinearProblem(a, L, bcs=[bc])
u_sol = problem.solve()

# Post-processing: Compute stress and strain
eps = epsilon(u_sol)
sig = sigma(u_sol)

# Save results
from dolfinx.io import XDMFFile

with XDMFFile(MPI.COMM_WORLD, "displacement.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(u_sol)

# For strain and stress, you may need to project them to a suitable function space
W = fem.TensorFunctionSpace(mesh, ("CG", 1))
strain = fem.project(eps, W)
stress = fem.project(sig, W)

with XDMFFile(MPI.COMM_WORLD, "strain.xdmf", "w") as file:
    file.write_function(strain)

with XDMFFile(MPI.COMM_WORLD, "stress.xdmf", "w") as file:
    file.write_function(stress)

print("FEM calculation complete. Results saved.")