How to interpolate the strain u as a Function(u) over the cross section nodes of the cantilever beam of hexahedrons?

General
1

dolfinx 0.8.0…

beam_mock.py:


import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem

L = 120.0         # Beam length in inches
W = 4.0           # Cross-sectional width (and height here; for a square beam)
E = 2e6     # Young's modulus in psi
nu = 0.3          # Poisson's ratio
mu = E / (2 * (1 + nu)) 
lambda_ = (E * nu) / ((1 + nu) * (1 - 2 * nu))
rho = 30.0 / 12**3  # Density in lb/in^3 (converted from 30 lb/ft^3)
domain = mesh.create_box(MPI.COMM_WORLD,
                         [np.array([0, 0, 0]), np.array([L, W, W])],
                         [121, 5, 5],
                         cell_type=mesh.CellType.hexahedron)
with io.XDMFFile(domain.comm, "beam_mesh.xdmf", "w") as xdmf_file:
    xdmf_file.write_mesh(domain)
w = -0.50
w_load = w * L
num_nodes = domain.topology.index_map(domain.topology.dim).size_local
load_per_node = w_load / num_nodes
domain = mesh.create_box(MPI.COMM_WORLD,
                         [np.array([0, 0, 0]), np.array([L, W, W])],
                         [121, 5, 5],
                         cell_type=mesh.CellType.hexahedron)

V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
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 = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D,
                     fem.locate_dofs_topological(V, fdim, boundary_facets),
                     V)
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
ds = ufl.Measure("ds", domain=domain)
def epsilon(u):
    return ufl.sym(ufl.grad(u))
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
def distribute_load(load_per_node, V):
    f = fem.Function(V)
    f.x.array[:] = 0.0  # Initialize to zero
    for i in range(num_nodes):
        f.x.array[i*3 + 2] = load_per_node  # Apply load in z-direction, accounting for 3 components per node
    return f
f = distribute_load(load_per_node*4, V)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
problem = LinearProblem(a, L, bcs=[bc],
                        petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
np.savetxt('U.csv', uh.x.array, delimiter=',')  # Save only the known DOFs
I=4*4**3/12
max_M = w*120.0**2/2
d_max = (w*120.0**4)/(8*E*I)
print(f'd_max:{d_max}')
print(f'minima: {min(uh.x.array)}')
print(f'max moment:{max_M}')



d_max:-0.30375
minima: -0.27241218582494253
max moment:-3600.0

mesh_dolfinxintegrate.py:


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


# Define beam parameters
L = 120.0         # Beam length in inches
W = 4.0           # Cross-sectional width (and height)
E = 2000000.0     # Young's modulus in psi
nu = 0.3          # Poisson's ratio
mu = E / (2 * (1 + nu))
lambda_ = (E * nu) / ((1 + nu) * (1 - 2 * nu))
rho = 30.0 / 12**3  # Density in lb/in^3
domain = mesh.create_box(MPI.COMM_WORLD,
                         [np.array([0, 0, 0]), np.array([L, W, W])],
                         [121, 5, 5],
                         cell_type=mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
U_array = np.loadtxt('U.csv', delimiter=',')
u = fem.Function(V)
u.x.array[:] = U_array

def strain(u):
    return 0.5 * (ufl.grad(u) + ufl.grad(u).T)
epsilon = strain(u)
def stress(epsilon, mu, lambda_):
    return 2 * mu * epsilon + lambda_ * ufl.tr(epsilon) * ufl.Identity(3)
sigma = stress(epsilon, mu, lambda_)

#       M_z = - ∫_A y σ_xx dA,
tol = 1e-8
def on_right(x):
    return np.isclose(x[0], L, atol=tol)
facet_indices = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, on_right)
num_facets = len(facet_indices)
tags = np.full(num_facets, 1, dtype=np.int32)
facet_tags = mesh.meshtags(domain, domain.topology.dim - 1, facet_indices, tags)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
x_coord = ufl.SpatialCoordinate(domain)
sigma_xx = sigma[0, 0]
Mz_expr = - (x_coord[1] - 2.0) * sigma_xx * ds(1)

Mz = fem.assemble_scalar(fem.form(Mz_expr))
if MPI.COMM_WORLD.rank == 0:
    print("Bending moment M_z at x = L is:", Mz)




Bending moment M_z at x = L is: 1.4755737098196857e-11

OK… So far M_z is much to small we should be -3500± inch/lbs of moment…I looked in ParaView to make sure that facet_indices are showing up and they are, so those entail the free end…

Why is M_z so miniscule? Can we get things to more accurately interpolate the situation?

2

you’re just measuring it along the wrong axis, it is x_coord[2] that you want if the load is acting in z direction. Also why are you applying a load like this by specifying values at the nodes? This does not corresponding to a correct constant distributed load because nodes may have a different number of incoming elements so, the equivalent nodal forces are not the same for all nodes. Just use a fem.Constant instead as in Implementation — FEniCSx tutorial

1 Like
6

I was able to get something remotely palatable to come out with a scale factor so far of only 175.00:

import numpy as np 
import meshio 
from includes.shape_hexahedron import get_unique_xs, cells_at_section, compute_hexahedron_volume, get_dofs_cell
from basix import ElementFamily, CellType, LagrangeVariant, create_element
from decimal import Decimal, getcontext
getcontext().prec = 150  # Set precision to 50 decimal places

def compute_strain_tensor(displacement_grad):
    strain = [[Decimal('0') for _ in range(3)] for _ in range(3)]
    for i in range(3):
        for j in range(3):
            strain[i][j] = Decimal('0.5') * (Decimal(str(displacement_grad[i][j])) + Decimal(str(displacement_grad[j][i])))
    return strain

def compute_stress_tensor(strain_tensor, lambda_, mu):
    stress = [[Decimal('0') for _ in range(3)] for _ in range(3)]
    trace = sum(Decimal(str(strain_tensor[i][i])) for i in range(3))
    for i in range(3):
        for j in range(3):
            stress[i][j] = lambda_ * trace * Decimal(str(int(i == j))) + Decimal('2') * mu * Decimal(str(strain_tensor[i][j]))
    return stress


def interpolate_displacement_gradient(gradients, epsilon):
    epsilon = epsilon.reshape(-1, 3)
    displacement_gradient = [[Decimal('0') for _ in range(3)] for _ in range(3)]
    for i, vertex_grad in enumerate(gradients):
        if i < len(epsilon):  # Ensure we don't go out of bounds
            for j in range(3):  # j is the component of displacement (u, v, or w)
                for k in range(3):  # k is the spatial direction (x, y, z)
                    if k < len(vertex_grad[0]):  # Check if k is within bounds of vertex_grad
                        displacement_gradient[k][j] += Decimal(str(epsilon[i, j])) * vertex_grad[0][k]
    return displacement_gradient


def create_gradient_table(mesh):
    # Create a Lagrange element of degree 1 (linear) for a hexahedron
    lagrange = create_element(ElementFamily.P, CellType.hexahedron, 1, LagrangeVariant.equispaced)
    gradients_table = {}
    for cell_index, cell_node_indices in enumerate(mesh.cells_dict['hexahedron']):
        cell_vertices = mesh.points[cell_node_indices]
        # Tabulate the derivatives at the vertices of the element
        tabulated = lagrange.tabulate(1, cell_vertices)
        # Extract gradients for each node of the element
        gradients = []
        for i in range(8):  # 8 nodes for hexahedron
            # For each node, we want the derivatives d/dx, d/dy, d/dz
            # Since tabulate gives [d/dx, d/dy, d/dz] for each node at each point, we reshape this:
            gradient = tabulated[1:, i, 0].reshape(3, 1)  # Grab the derivatives for node i
            # Convert the numpy array to a list of Decimal values
            gradients.append([[Decimal(str(g[0])) for g in gradient]])
        gradients_table[cell_index] = {'gradients': gradients}
    return gradients_table

# Main script
U, mesh = np.loadtxt('U.csv'), meshio.read('beam_mesh.xdmf')
U_decimal = [Decimal(str(u)) for u in U.flatten()]
hexahedron_indices = mesh.cells_dict['hexahedron']
unique_xs = get_unique_xs(mesh)
cells_incidental = cells_at_section(unique_xs[121], mesh)
E, nu, total_moment = Decimal('2e6'), Decimal('0.3'), np.zeros(3)
lambda_, mu = (E * nu) / ((Decimal('1') + nu) * (Decimal('1') - Decimal('2') * nu)), E / (Decimal('2') * (Decimal('1') + nu))
cell0_vertices = mesh.points[hexahedron_indices[0]]
print(f'cell0_vertices: {cell0_vertices}')
hexahedron_volume = compute_hexahedron_volume(cell0_vertices)
node_strain_tensors = {}
node_coordinates = {}
# Create the gradient table using basix
gradients_table = create_gradient_table(mesh)
for i, cell_index in enumerate(cells_incidental):
    cell_node_indices = hexahedron_indices[cell_index]
    if any(mesh.points[node][0] == unique_xs[121] for node in cell_node_indices):
        gradients = gradients_table[cell_index]['gradients']
        cell_dofs = get_dofs_cell(cell_node_indices)
        epsilon = U[cell_dofs]
        # Use epsilon directly without converting to local frame
        displacement_gradient = interpolate_displacement_gradient(gradients, epsilon)
        strain_tensor = compute_strain_tensor(displacement_gradient)
        for node in cell_node_indices:
            if node not in node_strain_tensors:
                node_strain_tensors[node] = []
            node_strain_tensors[node].append(strain_tensor)
            if node not in node_coordinates:
                node_coordinates[node] = mesh.points[node]
        stress_tensor = compute_stress_tensor(strain_tensor, lambda_, mu)
    else:
        pass
        #print(f"Skipping element {i} as no node has x = 120.0")
M = Decimal('0')
for node, strains in node_strain_tensors.items():
    for tensor in strains:
        c = Decimal(str(node_coordinates[node][2]))
        M += tensor[0][0] * c 
print(f'M:{M}')


# So this would need to be scaled up about 175.00x to meet the theoretical
# of -3600.00 inch/lbs
M:-20.8718885973339024957026409033302010804974220783706880 inch/lbs

@ x == 120.0 (the free end of the cantilever.)

Ok… I guess I realized that at this density and precision level that maybe the displacements would have to be scaled way up to accurately represent the scenario… Is that it? Is there more to it that I still don’t see it?

7

Please note that in all your posts, you keep shifting the issue/focus, and rarely address the questions that people ask you. For instance

It you are using this expression anywhere but at the end of the beam, it will give you zero. This is because ds

Is integrating over exterior facets, ie those that are connected to a single cell.

You also did not address

Finally, if you want to measure the stress as an interior cross section, keep in mind that sigma_xx will be discontinuous across the interior cell interfaces, and you should consistently pick a side. You can for instance do this with a custom integration domain, as shown in Add one-sided integration example as example of custom integration · Issue #158 · jorgensd/dolfinx-tutorial · GitHub

1 Like
8 

f = fem.Constant(domain, default_scalar_type((0, 0, -0.030)))

can be used to match-in to theoretical displacements. I had addressed it but it got deleted somewhere…


M_expr = - (x_coord[2]) * sigma_xx * ds(1)

Ok… I guess the above is what was meant… which also isn’t there through my struggle I had through all this…

Ok…


ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)

this above is exterior factets only… That is actually news to me at this point…

So you suggest:

so at this point I haven’t had an opportunity to look that over however the first thing that comes to mind is that if only the exterior facets are recovered there is there a way to recover the interior facets that are needed?

So far the reason for the tangents I had are that I was not yet aware that all of these things were going on…

9

Yes, that is the point of the reference. There you can use ‘ds’ for any facet, as it is specified by the cell index and local facet index.
In the example, it finds all cells one one side of an interface (based on cell markers). If one instead consider the midpoint of each cell, you can go from a facet on the interior, to find the cell to the left or right of this facet.

1 Like
10

I found the Jacobian to correct for large strains, from continuum mechanics, seems to work Ok… so far when used in conjunction with the area of the free end there… It’s being called F, or a part thereof…Comes within 3-500 inch/lbs of the theoretical…