Structural analysis of cantilever beam in multi-layer complexes

General
2

I reworked the code based on this article

> Blockquote

# 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)

T = functionspace(mesh, ("DG", 0, (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_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
#=======================

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


def left_boundary(x):
    return np.isclose(x[0], 0, atol=1e-6)



left_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, left_boundary)


mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)

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)

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.0001*9.81)))

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)

Now the displacement result is still not right，please.