Thank you for your reply, I’m sorry I didn’t articulate my question correctly. I created a model and divided the physics groups by using gmsh. Now the code works fine, but I’m not sure I’m correctly dispensing five different layers of material. Allow me to put the full code.
Here’s my code to create a model and divide the physical groups by using gmsh:
import gmsh
gmsh.initialize()
gmsh.model.add("CompositePlate")
# Define the overall box
dx, dy = 10., 2.
layer_heights = [0.2, 0.02, 0.2, 0.02, 0.2] #The thickness of each layer
dz = sum(layer_heights) # total thickness
layer_boxes = []
z_offset = 0
for i, thickness in enumerate(layer_heights):
box = gmsh.model.occ.addBox(0, 0, z_offset, dx, dy, thickness)
layer_boxes.append((3, box))
z_offset += thickness
gmsh.model.occ.fragment(layer_boxes, layer_boxes)
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(3)
for i, (dim, tag) in enumerate(volumes):
gmsh.model.addPhysicalGroup(dim, [tag], i + 1)
gmsh.model.setPhysicalName(dim, i + 1, f"Layer_{i+1}")
# Get All Edges (1D Segments)
entities = gmsh.model.getEntities(1)
edges_x = [] # X-direction
edges_y = [] # Y-direction
edges_z = [] # Z-direction
for dim, tag in entities:
xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(dim, tag)
if abs(xmax - xmin) > 1e-6 and abs(ymax - ymin) < 1e-6 and abs(zmax - zmin) < 1e-6:
edges_x.append(tag)
elif abs(xmax - xmin) < 1e-6 and abs(ymax - ymin) > 1e-6 and abs(zmax - zmin) < 1e-6:
edges_y.append(tag)
elif abs(xmax - xmin) < 1e-6 and abs(ymax - ymin) < 1e-6 and abs(zmax - zmin) > 1e-6:
edges_z.append(tag)
# Sets the number of cells in each direction
num_x = 50
num_y = 10
num_z_per_layer = 2 # 2 units in the Z direction of each floor
# Transfinite mesh
for edge in edges_x:
gmsh.model.mesh.setTransfiniteCurve(edge, num_x + 1)
for edge in edges_y:
gmsh.model.mesh.setTransfiniteCurve(edge, num_y + 1)
for edge in edges_z:
gmsh.model.mesh.setTransfiniteCurve(edge, num_z_per_layer + 1)
# Sets the faces to a structured quadrilateral mesh
surfaces = gmsh.model.getEntities(2) # Get All Faces
for dim, tag in surfaces:
gmsh.model.mesh.setTransfiniteSurface(tag)
gmsh.model.mesh.setRecombine(dim, tag)
# Set the volume to a structured hexahedral mesh
volumes = gmsh.model.getEntities(3) # Get All Volumes
for dim, tag in volumes:
gmsh.model.mesh.setTransfiniteVolume(tag)
# generate mesh
gmsh.model.mesh.generate(3)
gmsh.fltk.run()
# Export mesh file
mesh_file_name = "generated_mesh_1.msh"
gmsh.write(mesh_file_name)
print(f"Mesh successfully written to {mesh_file_name}")
gmsh.finalize()
Here’s the full code for solving with fenicsx:
import json
import numpy as np
from numpy.linalg import inv
import mpi4py
from dolfinx.mesh import locate_entities_boundary
from dolfinx import mesh, fem, plot, io, default_scalar_type
import ufl
import pyvista as pv
import dolfinx.fem.petsc
# Read mesh file
mesh, cell_tag, facet_tag = dolfinx.io.gmshio.read_from_msh(msh_path, mpi4py.MPI.COMM_WORLD, 0, gdim=3)
# Read the material properties JSON file
with open(json_path, 'r') as f:
materials = json.load(f)
all_layers = materials
# Rotation matrix
def rotation_matrix(C, theta):
# x_3 roate
theta = (np.pi * theta)
c = np.cos(theta)
s = np.sin(theta)
R = np.zeros((6, 6), dtype=np.float64)
R[0, 0] = c * c
R[0, 1] = s * s
R[0, 5] = 2 * c * s
R[1, 0] = s * s
R[1, 1] = c * c
R[1, 5] = -2 * c * s
R[2, 2] = 1
R[3, 3] = c
R[3, 4] = s
R[4, 3] = -s
R[4, 4] = c
R[5, 0] = -c * s
R[5, 1] = c * s
R[5, 5] = (c * c) - (s * s)
# Perform the rotation of the stiffness matrix C
Chat = np.matmul(R, np.matmul(C, np.transpose(R)))
return Chat
# Define the stiffness matrix of anisotropic materials
def get_elasticity_matrix(layer):
if layer["material_type"] == "Composite Layer (Anisotropic)":
E1, E2, E3 = layer["E1"], layer["E2"], layer["E3"]
nu21, nu31, nu32 = layer["nu21"], layer["nu31"], layer["nu32"]
G12, G13, G23 = layer["G12"], layer["G13"], layer["G23"]
fiber_angle = layer["fiber_angle"]
layer["fiber_angle"] = np.clip(perturbed_angle, 0, 1)
S = np.zeros((6, 6))
S[0, 0] = 1 / E1
S[0, 1] = -nu21 / E2
S[0, 2] = -nu31 / E3
S[1, 0] = S[0, 1]
S[1, 1] = 1 / E2
S[1, 2] = -nu32 / E3
S[2, 0] = S[0, 2]
S[2, 1] = S[1, 2]
S[2, 2] = 1 / E3
S[3, 3] = 1 / G23
S[4, 4] = 1 / G13
S[5, 5] = 1 / G12
# The Stiffness Matrix C is calculated as the inverse of the conformance matrix
C = inv(S)
elasticity_matrix_rotated = rotation_matrix(C, perturbed_angle)
return elasticity_matrix_rotated
else:
E, nu = layer["E"], layer["nu"]
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
elasticity_matrix = np.array([
[lmbda + 2 * mu, lmbda, lmbda, 0, 0, 0],
[lmbda, lmbda + 2 * mu, lmbda, 0, 0, 0],
[lmbda, lmbda, lmbda + 2 * mu, 0, 0, 0],
[0, 0, 0, mu, 0, 0],
[0, 0, 0, 0, mu, 0],
[0, 0, 0, 0, 0, mu]
])
return elasticity_matrix
# Obtain the stiffness matrix for all material layers
layer_materials = [get_elasticity_matrix(layer) for layer in all_layers]
#=======================
# assemble siffness tensor
#=======================
material_tags = np.unique(cell_tag.values)
# Define function space for different material
V_mat = functionspace(mesh, ("DG", 0, (6, 6)))
stiffness = dolfinx.fem.Function(V_mat)
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_tag.find(tag)
if tag <= len(layer_materials):
C = layer_materials[tag-1]
stiffness.interpolate(lambda x: tensor(x, C.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
#=======================
# Define Function Space
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))
# 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
# Define the judgment function of the left end face
def left_boundary(x):
return np.isclose(x[0], 0, atol=1e-6)
# Find Left Face
left_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, left_boundary)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
# Apply Fixed Boundary Condition (Left End Face)
left_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, left_facets)
u_bc = np.zeros(mesh.geometry.dim, dtype=np.float64)
bc_left = dolfinx.fem.dirichletbc(u_bc, left_dofs, V)
# Apply gravity
T = fem.Constant(mesh, default_scalar_type((0, 0, 0)))
ds = ufl.Measure("ds", domain=mesh)
f = fem.Constant(mesh, default_scalar_type((0, 0, -0.000981)))
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
# Solve
# Linear Problem
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc_left],
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
displacement = uh.x.array.reshape(-1, 3)
# Outputs the displacement of each layer
print("Layer\tMax Displacement\t\tMin Displacement\t\tMean Displacement")
for i in range(len(layer_materials)):
layer_cells = cell_tag.indices[cell_tag.values == (i + 1)]
layer_nodes = np.unique(np.concatenate(
[mesh.topology.connectivity(mesh.topology.dim, 0).links(c) for c in layer_cells]
))
layer_displacement = displacement[layer_nodes]
print(
f"{i + 1}\t"
f"{np.max(np.linalg.norm(layer_displacement, axis=1)):.6f}\t"
f"{np.min(np.linalg.norm(layer_displacement, axis=1)):.6f}\t"
f"{np.mean(np.linalg.norm(layer_displacement, axis=1)):.6f}"
)
I’m not sure if my question is posted here for consultation, but thanks again for your reply.
