How to set up beam kinematics for a squaure cross section of a solid beam of hexahedral elements?

General
1

https://bleyerj.github.io/comet-fenicsx/tours/beams/beams_3D/beams_3D.html

I have a modifed code:


import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem
import ufl
import basix 
import dolfinx.fem.petsc

comm = MPI.COMM_WORLD

# -------------------------------------------------------------------
# 1. Mesh and Domain Setup
# -------------------------------------------------------------------
# Define the dimensions of the box domain (beam)
Lx, Ly, Lz = 4.0, 8.0, 120.0  # Length, width, height
nx, ny, nz = 2, 2, 4  # Reduced mesh resolution


beam_width = 4.0
beam_height = 8.0
# Create a hexahedral mesh
domain = mesh.create_box(comm, [[0.0, 0.0, 0.0], [Lx, Ly, Lz]],
                         [nx, ny, nz], mesh.CellType.hexahedron)
gdim = domain.geometry.dim
tdim = domain.topology.dim

print(f"Geometrical dimension = {gdim}")
print(f"Topological dimension = {tdim}")

# -------------------------------------------------------------------
# 2. Create a UFL Function for the Geometry (Coordinate Mapping)
# -------------------------------------------------------------------
# In dolfinx 0.8.0, the correct function to create a function space is fem.functionspace
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
x_func = fem.Function(V)
# Interpolate the identity mapping: x = X
x_func.interpolate(lambda X: X)

dx_dX = ufl.Jacobian(domain)[:, 2]
t = dx_dX / ufl.sqrt(ufl.inner(dx_dX, dx_dX))
# Compute the unit tangent vector along the z-direction
t = dx_dX / ufl.sqrt(ufl.inner(dx_dX, dx_dX))

ez = ufl.as_vector([0, 0, 1])
a1 = ufl.cross(t, ez)
a1 /= ufl.sqrt(ufl.dot(a1, a1))
a2 = ufl.cross(t, a1)
a2 /= ufl.sqrt(ufl.dot(a2, a2))

# -------------------------------------------------------------------
# 3. Define Material and Cross-sectional Properties (4" x 8")
# -------------------------------------------------------------------
# Define cross-sectional dimensions (inches)
a = 4.0  # smaller side (inches)
b = 8.0  # larger side (inches)

# Compute the torsion constant J (approximation for a solid rectangular section)
# Formula: J ≈ (a * b^3 / 3) * [1 - 0.63*(a/b) + 0.052*(a/b)^5]
J_value = (a * b**3 / 3.0) * (1 - 0.63 * (a / b) + 0.052 * (a / b)**5)
print(f"Torsion constant J = {J_value:.2f} in^4")

# Material properties (in psi)
E_value = 70e3    # Young's modulus (psi)
nu_value = 0.3    # Poisson's ratio

# Create dolfinx constants for E and nu
E = fem.Constant(domain, E_value)
nu = fem.Constant(domain, nu_value)

# Compute shear modulus G = E / [2(1+nu)]
G_value = E_value / (2 * (1 + nu_value))
G = fem.Constant(domain, G_value)
print(f"Shear modulus G = {G.value:.2f} psi")

# Compute torsional stiffness: GJ = G * J
G = G_value * J_value
GJ = G_value * J_value

rho = fem.Constant(domain, 2.7e-3)
g = fem.Constant(domain, 9.81)

S = 4*8
ES = E*S
kappa = fem.Constant(domain, 5./6.)
GS1 = kappa*G*S
GS2 = kappa*G*S
EI1 = E*beam_width*beam_height**3/12
EI2 = E*beam_width**3*beam_height/12

Ue = basix.ufl.element("P", domain.basix_cell(), 1, shape=(gdim,))
W = fem.functionspace(domain, basix.ufl.mixed_element([Ue, Ue]))


u_, theta_ = ufl.TestFunctions(W)
du, dtheta = ufl.TrialFunctions(W)

def tgrad(u):
    return ufl.dot(ufl.grad(u), t)
def generalized_strains(u, theta):
    return ufl.as_vector([ufl.dot(tgrad(u), t),
                      ufl.dot(tgrad(u), a1)-ufl.dot(theta, a2),
                      ufl.dot(tgrad(u), a2)+ufl.dot(theta, a1),
                      ufl.dot(tgrad(theta), t),
                      ufl.dot(tgrad(theta), a1),
                      ufl.dot(tgrad(theta), a2)])
def generalized_stresses(u, theta):
    return ufl.dot(ufl.diag(ufl.as_vector([ES, GS1, GS2, GJ, EI1, EI2])), generalized_strains(u, theta))

Sig = generalized_stresses(du, dtheta)
Eps_ =  generalized_strains(u_, theta_)

dx_shear = ufl.dx(scheme="default",metadata={"quadrature_scheme":"default", "quadrature_degree": 0})

k_form = sum([Sig[i]*Eps_[i]*ufl.dx for i in [0, 3, 4, 5]]) + (Sig[1]*Eps_[1]+Sig[2]*Eps_[2])*dx_shear
l_form = -rho*S*g*u_[2]*ufl.dx


def left_end(x):
    return np.isclose(x[2], 0.0)

def right_end(x):
    return np.isclose(x[2], Lz)

Vu, _ = W.sub(0).collapse()
Vt, _ = W.sub(1).collapse()

# Locate degrees of freedom at the left end (z=0) and right end (z=Lz)
u_dofs_left = fem.locate_dofs_geometrical((W.sub(0), Vu), left_end)
theta_dofs_left = fem.locate_dofs_geometrical((W.sub(1), Vt), left_end)

u_dofs_right = fem.locate_dofs_geometrical((W.sub(0), Vu), right_end)
theta_dofs_right = fem.locate_dofs_geometrical((W.sub(1), Vt), right_end)

# Functions for the Dirichlet boundary conditions
u0 = fem.Function(Vu)
theta0 = fem.Function(Vt)

# Apply boundary conditions to both ends
bcs = [fem.dirichletbc(u0, u_dofs_left, W.sub(0)),
       fem.dirichletbc(theta0, theta_dofs_left, W.sub(1)),
       fem.dirichletbc(u0, u_dofs_right, W.sub(0)),
       fem.dirichletbc(theta0, theta_dofs_right, W.sub(1))]

w = fem.Function(W, name="Generalized_displacement")

problem = fem.petsc.LinearProblem(
    k_form, l_form, u=w, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
problem.solve()

print('done...')

So there is a tangent vector to use… OK… I tried to modify that to match the condition now as opposed to what is in the original example.

I made these modifications in an attempt to meet the current conditions:


# Compute shear modulus G = E / [2(1+nu)]
G_value = E_value / (2 * (1 + nu_value))
G = fem.Constant(domain, G_value)
print(f"Shear modulus G = {G.value:.2f} psi")
# Compute torsional stiffness: GJ = G * J
G = G_value * J_value
GJ = G_value * J_value
rho = fem.Constant(domain, 2.7e-3)
g = fem.Constant(domain, 9.81)
S = 4*8
ES = E*S
kappa = fem.Constant(domain, 5./6.)
GS1 = kappa*G*S
GS2 = kappa*G*S
EI1 = E*beam_width*beam_height**3/12
EI2 = E*beam_width**3*beam_height/12

So far I have this traceback to deal with:


Geometrical dimension = 3
Topological dimension = 3
Torsion constant J = 468.74 in^4
Shear modulus G = 769230.77 psi
Traceback (most recent call last):
  File "/home/prusso/dolfinx-rotations/example003.py", line 142, in <module>
    problem = fem.petsc.LinearProblem(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 762, in __init__
    self._a = _create_form(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 249, in form
    return _create_form(form)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 244, in _create_form
    return _form(form)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 186, in _form
    ufcx_form, module, code = jit.ffcx_jit(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/jit.py", line 51, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/jit.py", line 201, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 276, in compile_forms
    raise e
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 256, in compile_forms
    impl = _compile_objects(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 383, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/compiler.py", line 113, in compile_ufl_objects
    ir = compute_ir(analysis, _object_names, _prefix, options, visualise)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/ir/representation.py", line 224, in compute_ir
    irs = [
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/ir/representation.py", line 225, in <listcomp>
    _compute_integral_ir(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/ir/representation.py", line 588, in _compute_integral_ir
    integral_ir = compute_integral_ir(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/ir/integral.py", line 85, in compute_integral_ir
    mt_table_reference = build_optimized_tables(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/ir/elementtables.py", line 427, in build_optimized_tables
    tabletype = analyse_table_type(tbl)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/ir/elementtables.py", line 576, in analyse_table_type
    elif is_quadrature_table(table, rtol=rtol, atol=atol):
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/ir/elementtables.py", line 538, in is_quadrature_table
    Id = np.eye(num_points)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/numpy/lib/_twodim_base_impl.py", line 219, in eye
    m = zeros((N, M), dtype=dtype, order=order, device=device)
numpy._core._exceptions._ArrayMemoryError: Unable to allocate 1.40 TiB for an array with shape (438976, 438976) and data type float64
Exception ignored in: <function LinearProblem.__del__ at 0x74d225bd05e0>
Traceback (most recent call last):
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 805, in __del__
    self._solver.destroy()
AttributeError: 'LinearProblem' object has no attribute '_solver'

It seems to state numpy._core._exceptions._ArrayMemoryError: Unable to allocate 1.40 TiB for an array with shape (438976, 438976) and data type float64… OK… I tried to reduce the resolution of the box mesh as much as I could however I couldn’t get the traceback to go away…

The traceback seems rather errouneous in nature for a mesh this size and resolution…

Did these modifications I made to the tangent vector and B matrix adequetly match the conditions that are now in place?

*Why is the traceback appearing? *

How can I get the traceback to go away?

2

Why do you use a hexahedral mesh if you want to solve a beam problem? Beam theory is formulated on the 1d domain corresponding to the beam median axis. If you want to solve the deformation of a beam-like geometry using solid 3d elements you don’t need the rotation degrees of freedom theta. Just use the standard 3d equations of elasticity using the displacement field.

1 Like
3

The use of the hexahedral elements has to do with the generation of a Stress map through the whole solid body. Later on, what I could be looking to do is mock-up solid of the rough dimensions with multiple materials. OK, I accept good advice here.

There is a publication out there how to set up Hermite hexahedron to allow for the specific entry of 6 degrees of freedom per node that does allow for the specific entry of (fx,fy,fz,mx,my,mz) and the solving with an assembled K matrix…

Doing so for the square beam problem, was a can of worms that I know of. The scenario leads to a singularity, or a situation where the K need a large regularization parameter placed on the diagonal of the K for a well distributed stress map to come out, and the regularization parameter has to be chosen perfectly for the specific mock otherwise the stress distribution comes out no good. Although so far as I know if all of these quirks and nuances with difficulties are dealt with for that Hermite method (B matrix) the stress maps are accurate and very beautiful…

I’ll bet that is because I just don’t really know yet everything and I’m nudging the problem wrongly somewhere in my ignorance. The displacements are so far very small in compared to theoretical values of such a Hermite formulation after things are regularized however if the regularization parameter and the scaling of the shape functions Hn1-16 are perfect the stress does come out within about 20% of what things should be…

I have a great appreciation that you would write me back… So as I have been looking for a better way to do it and it has been a prolonged turmoil for me, to get out of the dark. OK… I really like:

https://jsdokken.com/dolfinx-tutorial/chapter2/linearelasticity.html

so I guess that using Jørgen’s is what you suggest. OK with me… So in order to get all this right I will need to integrate the surface that represents the cross section over the whole domain to get the moment, which will be important to match up mocks to the building code…

Sorry, I haven’t meant to be a bother… it’s just that getting all this as well done as I can is a thing I have been meaning to do for a very long time. I didn’t know how much I was taking on when I started…Thanks so much! Both you and Jørgen have been a super help, so without the acceptance of guidance here I would be very in the dark still.