Thank you so much!
For everyone interested, here is the complete code:
from pymaterial.materials import IsotropicMaterial, TransverselyIsotropicMaterial, OrthotropicMaterial
import numpy as np
import numpy as np
import gmsh
from mpi4py import MPI
mat0 = {
"thickness": 0.3,
"stiff_tens": IsotropicMaterial(Em=2e11,nu=0.3, density=7850).get_stiffness() # steel from ANSYS
}
mat1 = {
"thickness": 0.3,
"stiff_tens": OrthotropicMaterial(E_x=4.5e10,
E_y=1e10,
E_z=1e10,
nu_xy=0.3,
nu_yz=0.4,
nu_xz=0.3,
G_xy=5e9,
G_yz=3.8462e9,
G_xz=5e9,
density=2000
).get_stiffness() # from ANSYS Epoxy E-Glass UD
}
mat2 = {
"thickness": 0.3,
"stiff_tens": OrthotropicMaterial(E_x=2.3e11,
E_y=2.3e10,
E_z=2.3e10,
nu_xy=0.2,
nu_yz=0.4,
nu_xz=0.2,
G_xy=9e9,
G_yz=8.2143e9,
G_xz=9e9,
density=1800
).get_stiffness() # from ANSYS Carbon Fiber (230 GPa)
}
stack = [
mat0,
mat1,
mat2
]
width, length = 1, 1
thickness = 0
gdim = 3 # Geometric dimension of the mesh
model_rank = 0
mesh_comm = MPI.COMM_WORLD
#=================
# generate mesh
#=================
gmsh.initialize()
model = gmsh.model()
name = "Battery Stack"
model.add(name)
model.setCurrent(name)
if mesh_comm.rank == model_rank:
boxes = []
coms = []
for i in range(len(stack)):
elem = stack[i]
elem_thickness = elem["thickness"]
start_height = thickness
end_height = thickness + elem_thickness
boxes.append(model.occ.addBox(0, 0, start_height, width, length, elem_thickness))
thickness = end_height
coms.append([width/2, length/2, start_height + elem_thickness/2])
model.occ.synchronize()
for i in range(len(boxes)-1):
results = model.occ.fragment([(3,boxes[i])], [(3,boxes[i+1])])
print(results)
model.occ.synchronize()
for dim_,id_ in model.getEntities(3):
com = gmsh.model.occ.getCenterOfMass(dim_, id_)
# print(com)
for i in range(len(coms)):
if not np.allclose(com, coms[i]):
continue
tag = i
model.addPhysicalGroup(3, [id_], tag=tag)
model.setPhysicalName(3, tag, stack[i]["name"])
break
model.mesh.generate(dim=3)
# gmsh.model.mesh.optimize("Netgen")
# Create DOLFINx distributed mesh
from dolfinx.io.gmshio import model_to_mesh
mesh, cell_markers, facet_markers = model_to_mesh(model, mesh_comm, model_rank)
gmsh.finalize()
#=======================
# assemble siffness tensor
#=======================
material_tags = np.unique(cell_markers.values)
T = dolfinx.fem.TensorFunctionSpace(mesh, ("DG", 0), shape=(6, 6))
stiffness = dolfinx.fem.Function(T)
def tensor(x, tens):
tens_reshaped = tens.reshape((6 * 6, 1))
return np.broadcast_to(tens_reshaped, (6 * 6, x.shape[1]))
for tag in material_tags:
domain = cell_markers.find(tag)
if tag < len(stack):
material = stack[tag]
stiffness.interpolate(lambda x: tensor(x, material["stiff_tens"].astype(np.float32)), domain)
#=======================
# helpers
#=======================
def epsilon(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def strain2voigt(e, dim: int = 3):
"""e is a 2nd-order tensor, returns its Voigt vectorial representation"""
if dim == 3:
return ufl.as_vector([e[0,0],e[1,1],e[2,2],2*e[1,2], 2*e[0,2],2*e[0,1]])
else:
return ufl.as_vector([e[0,0],e[1,1],2*e[0,1]])
def voigt2stress(s, dim: int = 3):
"""
s is a stress-like vector (no 2 factor on last component)
returns its tensorial representation
"""
if dim == 3:
return ufl.as_tensor([
[s[0], s[5], s[4]],
[s[5], s[1], s[3]],
[s[4], s[3], s[2]]
])
else:
return ufl.as_tensor([
[s[0], s[2]],
[s[2], s[1]]
])
def sigma(u, stiff_tens):
return voigt2stress(ufl.dot(stiff_tens, strain2voigt(epsilon(u))))
#=======================
# problem
#=======================
from dolfinx.fem import FunctionSpace
import ufl
element = ufl.VectorElement('CG', mesh.ufl_cell(), 2)# CG-Lagrange element, 2-degree of freedom (2nd order element)
V = FunctionSpace(mesh, element)
# Define variational problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(
ufl.dot(ufl.as_matrix(stiffness), strain2voigt(epsilon(u))),
strain2voigt(epsilon(v))
) * ufl.dx
#=======================
# loading
#=======================
from dolfinx.fem import Constant
from petsc4py.PETSc import ScalarType
import ufl
body_force = Constant(mesh, ScalarType((0, 0, 0))) # force on body
def surface_force(x):
return np.isclose(x[1], length)
fdim = mesh.topology.dim - 1
traction_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, surface_force)
traction = Constant(mesh, ScalarType((0, 10/((0.3*3)*width), 0))) # force on surface
ds = ufl.Measure("ds", domain=mesh)
L = ufl.inner(body_force, v) * ufl.dx + ufl.inner(traction, v) * ds
#=======================
# boundary conditions
#=======================
def clamped_boundary(x):
return np.isclose(x[1], 0)
fdim = mesh.topology.dim - 1
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, clamped_boundary)
u_D = np.array([0,0,0], dtype=ScalarType)
bc = dolfinx.fem.dirichletbc(u_D, dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets), V)
#=======================
# solve
#=======================
from dolfinx.fem.petsc import LinearProblem
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
#=======================
# export
#=======================
from dolfinx import io
from pathlib import Path
with io.XDMFFile(mesh.comm, Path("/home/fenics/shared/") / "output.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(cell_markers)
uh.name = "Deformation"
xdmf.write_function(uh)