I was able to make some improvements to the script so far. There is some scalar value of output for the Mz of -473±, however that seems to me much too small so far to be within an acceptable tolerance to what solid mechanics may predict.

```
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
# Geometry and material properties
L = 120.0
W = 4.0
H = 8.0
rho = 0.48
delta = W / L
gamma = 0.4 * delta**2
beta = 1.25
E = 2000000 # Modulus of elasticity
nu = 0.3 # Poisson's ratio
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
g = gamma
# Create the mesh
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, H])],
[20, 60, 60], cell_type=mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.topology.dim, )))
# Define boundary conditions
def clamped_boundary(x):
return np.isclose(x[0], 0)
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = fem.Constant(domain, np.array([0, 0, 0], dtype=default_scalar_type))
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)
def free_end(x):
return np.isclose(x[0], L)
# Find facets on the boundary where x = L (free end)
fdim = domain.topology.dim - 1
# Ensure that the facets at the free end (x = L) are correctly marked
free_end_facets = mesh.locate_entities_boundary(domain, fdim, free_end)
facet_markers = mesh.meshtags(domain, fdim, free_end_facets, np.full(len(free_end_facets), 2, dtype=np.int32))
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)
# Define the elasticity problem
def epsilon(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
fx = 0.1
fy = -rho * g * 0.25
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, fy, 0.0)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# Calculate bending moment Mz
x = ufl.SpatialCoordinate(domain)
Mz = fem.assemble_scalar(fem.form(-x[1] * sigma(uh)[0, 0] * ds()))
print(f"Bending moment Mz = {Mz:.6f}")
# Bending moment Mz = -473.749420
```

The density of the mesh so far is approaching the limitations of what I can accomplish locally on my own system. I noticed that if I increase density much further I get a message output to the terminal that states, “killed”, (more information about the message would be great!).

I don’t quite understand it at this point why -473, that seems to be a large departure from the 2D sample which had looked close enough for all intensive purposes to what solid mechanics may predict…

I noticed so far that the in terms of the 2D samples that the virtual work method does seem to come closer to what solid mechanics may predict. At this point I don’t really have a great handle on the virtual work method and how to employ the use of it outside on the context of the 2D sample that is shown…

Anyone can post this thread, why the Mz that is output here -473± seems to be roughly half of what is predicted in the 2D sample I provided a link to here on this thread is stating?