Why is 'fem.assemble_scalar(fem.form(ufl.cross(ufl_vec, sigma(u)*n)[1]*ds(ds_section)))' throwing "ValueError: Index out of bounds."

General
1

** **dolfinx 0.8.0...** **

From an earlier posting:

There is a reference to:

M=∫Sx×σn dS


n = ufl.FacetNormal(domain)
My = fem.assemble_scalar(fem.form(ufl.cross(x, sigma(u)*n)[1]*ds(surf_idx)))

Here is the full code I have:


import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import fem, mesh, io
import dolfinx.fem.petsc
from dolfinx.fem.petsc import LinearProblem
import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation


# Define the beam properties
L = 120.0  # Length in inches
H = 8.0    # Height in inches
W = 4.0    # Width in inches
I = (W * H**3) / 12  # Moment of inertia in inches^4
w = -0.5

# Mesh properties
Nx = 121  # Finer mesh along the length
Ny = 9  # Finer mesh along the height

print('NX:', Nx)
print('Ny:', Ny)

# Create the mesh
domain = mesh.create_rectangle(
    MPI.COMM_WORLD,
    [(0.0, -H / 2), (L, H / 2)],
    (Nx, Ny),
    diagonal=mesh.DiagonalType.crossed,
)
# Material properties
E = fem.Constant(domain, 2000000.0)
nu = fem.Constant(domain, 0.3)
mu = (E / (2 * (1 + nu)))
lmbda = (E * nu) / ((1 + nu) * (1 - 2 * nu))
points = domain.geometry.x
unique_x_values = np.unique(points[:, 0])  


# Define strain and stress
def eps(v):
    return ufl.sym(ufl.grad(v))

def sigma(v):
    return lmbda * ufl.tr(eps(v)) * ufl.Identity(2) + 2.0 * mu * eps(v)

# Applied load
fx = 0.0
fy = (w * L) / (Nx * Ny)
fy_analytic = 0.25
f = fem.Constant(domain, (fx, -fy))  # Uniformly distributed load over the length

# Define function space
V = fem.functionspace(domain, ("P", 2, (2,)))

# Define variational problem
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
a = ufl.inner(sigma(du), eps(u_)) * ufl.dx
l = ufl.inner(f, u_) * ufl.dx

# Boundary conditions
def left_end(x):
    return np.isclose(x[0], 0.0)


# Define the free-end function
def free_end(x):
    return np.isclose(x[0], L)

# fdim is the codimension (domain.topology.dim - 1)
fdim = domain.topology.dim - 1

# Locate the free-end facets
facets_free_end = mesh.locate_entities_boundary(domain, fdim, free_end)
# Create an array of markers for these facets. Here, we use the marker value 2.
markers_free_end = np.full(facets_free_end.shape, 2, dtype=np.int32)
# Create the meshtags for the free-end facets
facettags_free_end = mesh.meshtags(domain, fdim, facets_free_end, markers_free_end)
# Now define the integration measure using these meshtags
ds_free_end = ufl.Measure("ds", domain=domain, subdomain_data=facettags_free_end)

# Locate facets on the left (fixed end)
facets_left_end = mesh.locate_entities_boundary(domain, fdim, left_end)
markers_left_end = np.full_like(facets_left_end, 1)
facettags_left_end = mesh.meshtags(domain, fdim, facets_left_end, markers_left_end)
# Measure for boundary integration
ds_leftend = ufl.Measure("ds", domain=domain, subdomain_data=facettags_left_end)

# Apply Dirichlet BCs (fixing the left end)
left_dofs = fem.locate_dofs_geometrical(V, left_end)
u_bc = fem.Function(V)
bcs = [fem.dirichletbc(u_bc, left_dofs)]

# Solve the problem
u = fem.Function(V, name="Displacement")
problem = LinearProblem(a, l, bcs, u=u)
problem.solve()

x = ufl.SpatialCoordinate(domain)
n = ufl.FacetNormal(domain)
ds = ufl.Measure("ds", domain=domain)

#My = fem.assemble_scalar(fem.form(ufl.cross(x, sigma(u)*n)[1]*ds(ds_free)))
#My = fem.assemble_scalar(fem.form(ufl.cross(x, sigma(u)*n)[1]*ds(surf_idx)))

# Define points along the beam length where you want to evaluate M(x)

bending_moments = []


M_max = (w*L**2)/(2)
print(f'M_max: {M_max}')

for x_pos in unique_x_values:
    # Define a broad tolerance
    delta = 2 * L / (Nx - 1)  # Twice the element size along x-direction

    def cells_at_xpos(x):
        return np.isclose(x[0], x_pos, atol=delta)

    cells = mesh.locate_entities(domain, domain.topology.dim, cells_at_xpos)
    print(cells)

    if len(cells) == 0:
        print(f"No cells found at x = {x_pos}")
        bending_moments.append(0)
        continue  # Skip to next x_pos

    # Create cell tags and measure
    cell_indices = np.full(len(cells), 4, dtype=np.int32)  # Use marker '4' for these cells
    cell_tags = mesh.meshtags(domain, domain.topology.dim, cells, cell_indices)
    ds_section = ufl.Measure("dx", domain=domain, subdomain_data=cell_tags)

    # Define spatial coordinates and stress
    x = ufl.SpatialCoordinate(domain)
    n = ufl.FacetNormal(domain)
    sgn = sigma(u)*n
    ufl_vec = ufl.as_vector([x_pos,3.9])
    ufl_cross = ufl.cross(ufl_vec, sgn) # (3,)
    scalar_cross = ufl_cross[1]

    '''
    The following line will result in a ValueError: Index out of bounds.
    '''
    My = fem.assemble_scalar(fem.form(ufl.cross(ufl_vec, sigma(u)*n)[1]*ds(ds_section)))
    print(My*W)


# Plot bending moment along the beam length
import matplotlib.pyplot as plt

#plt.plot(x_positions, bending_moments)
#plt.xlabel('Position along beam length (inches)')
#plt.ylabel('Bending Moment (in-lb)')
#plt.title('Bending Moment Distribution')
#plt.grid(True)
#plt.show()

This is the problem line:


My = fem.assemble_scalar(fem.form(ufl.cross(ufl_vec, sigma(u)*n)[1]*ds(ds_section)))

Here is the resulting “ValueEror: Index out of bounds.” traceback:


fem.assemble_scalar(fem.form(ufl.cross(ufl_vec, sigma(u)*n)[1]*ds(ds_section)))
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  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 108, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"])  # type: ignore
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/analysis.py", line 94, in analyze_ufl_objects
    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/analysis.py", line 94, in <genexpr>
    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ffcx/analysis.py", line 180, in _analyze_form
    form_data = ufl.algorithms.compute_form_data(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/algorithms/compute_form_data.py", line 278, in compute_form_data
    form = preprocess_form(form, complex_mode)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/algorithms/compute_form_data.py", line 232, in preprocess_form
    form = apply_algebra_lowering(form)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/algorithms/apply_algebra_lowering.py", line 150, in apply_algebra_lowering
    return map_integrand_dags(LowerCompoundAlgebra(), expr)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/algorithms/map_integrands.py", line 86, in map_integrand_dags
    return map_integrands(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/algorithms/map_integrands.py", line 29, in map_integrands
    mapped_integrals = [
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/algorithms/map_integrands.py", line 30, in <listcomp>
    map_integrands(function, itg, only_integral_type) for itg in form.integrals()
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/algorithms/map_integrands.py", line 39, in map_integrands
    return itg.reconstruct(function(itg.integrand()))
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/algorithms/map_integrands.py", line 87, in <lambda>
    lambda expr: map_expr_dag(function, expr, compress), form, only_integral_type
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag
    (result,) = map_expr_dags(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/corealg/map_dag.py", line 103, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/algorithms/apply_algebra_lowering.py", line 60, in cross
    return as_vector((c(1, 2), c(2, 0), c(0, 1)))
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/algorithms/apply_algebra_lowering.py", line 58, in c
    return Product(a[i], b[j]) - Product(a[j], b[i])
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/exproperators.py", line 383, in _getitem
    all_indices, slice_indices, repeated_indices = create_slice_indices(
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/index_combination_utils.py", line 150, in create_slice_indices
    raise ValueError("Index out of bounds.")
ValueError: Index out of bounds.

So far I am looking for outputs on or near M_max: -3600.0, or less so…

Why is this happening? How to resolve this?

2

Probably because ufl.cross is only implemented for 3D vectors, ref: Form language — Unified Form Language (UFL) 2024.2.0 documentation

Note that this operator is only defined for vectors of length three.

1 Like
3 

    ufl_vec = ufl.as_vector([x_pos,3.9])
    ufl_cross = ufl.cross(ufl_vec, sgn) # The shape is: (3,)
    scalar_cross = ufl_cross[1]

OK… Here above, ufl.cross is used in conjuction with a two component ufl vector and the system will accept it…The shape is a 3 component vector, I guess… There is no way to print anything out, things are symbolic…

~~~
# I modified this in an attempt to retrieve edges, using that degree of freedom...
cells = mesh.locate_entities(domain, domain.topology.dim - 1, cells_at_xpos)
~~~
My = fem.assemble_scalar(fem.form(ufl.cross(ufl_vec, sigma(u)*n)[1]*ds(ds_section)))
~~~

So, ds_section should represent edges I guess, and ds should represent facets of the whole surface I guess… So there is a 3-component vector, where Jeremy says to use component [1] or the y-component… that 3-component vector is being multiplied by ds(ds_section)…

As far as ds(ds_section) what are the ramifications of doing this, or in other words what does this mean or represent? So far the actual value of ds(ds_section) is symbolic and there is no way here to print an actural value…

Even I was able to surmise so far that “*ds(ds_section)” in conjunction with that 3-component vector is an incompatible scenario…

While your reply so far is very helpful, I still don’t have really much to go with as far as how to fix that… although I tried to get a ufl.shape of the resultant ds(ds_section), things were not really a correct usage there or there is no specific shape to grab…

I also attempted a few times a conversion into 3D space of the whole transaction however I wasn’t yet able to accomplish such a thing…

Could I please get a bit more to run from here in terms of solving this? How is it possible to complete such a transaction in a compatible way?

4

What is unclear to me what you can extract from the [1] component of a cross-product of a 2D vector. It would at most be the third component of the cross product that would be interesting, ref: https://www.quora.com/Can-you-cross-product-2D-vectors.

What I would assume, is that you can implement the 2D cross product yourself:

   def cross_2D(x, y):
        return x[0] * y[1] - x[1] * y[0]
    sgn = sigma(u) * n
    ufl_vec = ufl.as_vector([x_pos, 3.9])
    My = fem.assemble_scalar(fem.form(cross_2D(ufl_vec, sigma(u) * n) * ds(ds_section)))

However, your code is rather long, and doesn’t get straight to the point, which is:
How to compute a cross product of 2D vectors with ufl.

1 Like