Shape mismatch building linear form for a mixed function space structural problem with 6 degrees of freedom per node

[prusso630]

This is the whole code:

import dolfinx
from dolfinx import mesh, fem, plot
from mpi4py import MPI
import basix
import numpy as np
from petsc4py import PETSc
import ufl
from dolfinx.fem.petsc import LinearProblem
from dolfinx.fem import dirichletbc, locate_dofs_topological
from ufl import Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, inner, dot
import pyvista as pv

# Material properties
E = 2000000.0  # Young's modulus in PSI (for structural wood)
nu = 0.299999  # Poisson's ratio (for structural wood)
L = 72  # Beam length in inches
W = 4.0  # Beam width in inches
w = 100.0  # Distributed load in lbs/inch

# Define the mesh with tetrahedral cells
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
                         [20, 6, 6], cell_type=mesh.CellType.tetrahedron)

# Define elements for displacements (3 DOF) and rotations (3 DOF)
elem_disp = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(3,))
elem_rot = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(3,))
elem_mixed = basix.ufl.mixed_element([elem_disp, elem_rot])

V = fem.functionspace(domain, elem_mixed)
VV = fem.functionspace(domain, elem_disp)

# Strain tensor as a symmetric gradient of the displacement field
def epsilon(u):
    return 0.5 * (ufl.grad(u) + ufl.grad(u).T)

# Stress tensor derived from Hooke's law
def sigma(u):
    lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))  # First Lamé parameter
    mu = E / (2 * (1 + nu))  # Shear modulus (second Lamé parameter)
    return lambda_ * ufl.tr(epsilon(u)) * ufl.Identity(3) + 2 * mu * epsilon(u)

# Define trial and test functions for displacements (u) and rotations (theta)
(u, theta) = TrialFunctions(V)
(v, tau) = TestFunctions(V)

# Weak form for the elasticity problem (displacement part)
a_disp = inner(sigma(u), epsilon(v)) * ufl.dx

# Rotational-related terms (rotations theta and tau)
a_rot = inner(theta, tau) * ufl.dx

# Coupling terms between displacement and rotation (optional)
a_coupling = inner(theta, v) * ufl.dx + inner(u, tau) * ufl.dx

# Final variational form: combining displacements, rotations, and coupling
a = a_disp + a_rot + a_coupling

# Get subspace of V for displacements
V0 = V.sub(0)
Q, _ = V0.collapse()
fdim = domain.topology.dim - 1

# Locate the entities where x=0 (left boundary) and x=L (right boundary)
def locate_left_boundary(x):
    return np.isclose(x[0], 0.0)

def locate_right_boundary(x):
    return np.isclose(x[0], L)

facets_left = mesh.locate_entities_boundary(domain, fdim, locate_left_boundary)
facets_right = mesh.locate_entities_boundary(domain, fdim, locate_right_boundary)

# Apply boundary condition at the left boundary where x = 0 (clamping)
dofs_left = locate_dofs_topological(V0, fdim, facets_left)
u_left = fem.Function(Q)
with u_left.vector.localForm() as loc:
    loc.set(0.0)  # Zero displacement in all directions
bc_left = dirichletbc(u_left, dofs_left)

# Apply boundary condition at the right boundary where x = L (clamping)
dofs_right = locate_dofs_topological(V0, fdim, facets_right)
u_right = fem.Function(Q)
with u_right.vector.localForm() as loc:
    loc.set(0.0)  # Zero displacement in all directions
bc_right = dirichletbc(u_right, dofs_right)

# Compute the distributed load per node and apply it
nodes = domain.geometry.x
num_nodes = len(nodes)

# Create force function f in the same mixed space V
f = fem.Function(V)

# Compute the distributed load per node and apply it
nodes = domain.geometry.x
num_nodes = len(nodes)

# Assign distributed force (downward in -y) and moment (around x-axis)
I = (W**4) / 12  # Moment of inertia for a rectangular cross-section

# Compute the number of degrees of freedom per node
dofs_per_node = 6

# Loop over all nodes to apply the distributed load and moment
for i, node in enumerate(nodes):
    fy = -w * L / num_nodes  # Distributed load down the -y axis
    y = node[1]  # y-coordinate of the node
    mx = -(w * abs(y) * L**2) / (2 * I)  # Moment around x-axis

    # Apply fy to the displacement component in y and mx to the rotational component around x
    f.vector.setValue(i * dofs_per_node + 1, fy)  # Displacement force in y
    f.vector.setValue(i * dofs_per_node + 3, mx)  # Rotational moment around x

# Define the right-hand side (RHS) for the variational problem
l = inner(f, v) * ufl.dx  # Adjust RHS to match the shape of the test function v

# Create the linear problem and solve it
#problem = LinearProblem(a, l, bcs=[bc_left, bc_right])
#uh = problem.solve()

# Save the solution to a CSV file
#np.savetxt('U.csv', uh.x.array)
#print('done...')

with traceback:

Traceback (most recent call last):
  File "/home/prusso/dolfinx-beam_w/compute_displacements2.py", line 115, in <module>
    l = inner(f, v) * ufl.dx  # Adjust RHS to match the shape of the test function v
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/operators.py", line 213, in inner
    return Inner(a, b)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/tensoralgebra.py", line 166, in __new__
    raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
ValueError: Shapes do not match: <Coefficient id=128873217154304> and <ListTensor id=128873219570112>

If I understand it, there is a shape problem building the linear form after using the force vector this way.

This is the desired algorithm so far:
(Set a distributed load, small cursive w, and a moment of force for each node):

# Loop over all nodes to apply the distributed load and moment
for i, node in enumerate(nodes):
    fy = -w * L / num_nodes  # Distributed load down the -y axis
    y = node[1]  # y-coordinate of the node
    mx = -(w * abs(y) * L**2) / (2 * I)  # Moment around x-axis

    # Apply fy to the displacement component in y and mx to the rotational component around x
    f.vector.setValue(i * dofs_per_node + 1, fy)  # Displacement force in y
    f.vector.setValue(i * dofs_per_node + 3, mx)  # Rotational moment around x

As an alternative to troubleshoot the problem I had also tried this way:

# Get subspace of V for displacements
V0 = V.sub(0)
Q, _ = V0.collapse()
fdim = domain.topology.dim - 1


# Compute the distributed load per node and apply it
nodes = domain.geometry.x
num_nodes = len(nodes)

f_disp = fem.Function(Q)

# Loop over all nodes to apply the distributed load (displacement force)
for i, node in enumerate(nodes):
    fy = -w * L / num_nodes  # Distributed load down the -y axis
    f_disp.vector.setValue(i * 3 + 1, fy)  # Displacement force in y

# Define the right-hand side (RHS) for the variational problem
l = inner(f_disp, v[0]) * ufl.dx  # Adjust RHS to match the shape of the test function v

which leads to a different traceback:

Traceback (most recent call last):
  File "/home/prusso/dolfinx-beam_w/compute_displacements3.py", line 99, in <module>
    l = inner(f_disp, v[0]) * ufl.dx  # Adjust RHS to match the shape of the test function v
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/operators.py", line 213, in inner
    return Inner(a, b)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/tensoralgebra.py", line 166, in __new__
    raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
ValueError: Shapes do not match: <Coefficient id=132880029505664> and <Indexed id=132880029444000>

What is going wrong here? How can I set up the desired algorithm to work as far as setting the forces and moments for this type of problem?

[prusso630]

Sorry about all this I was all thumbs at first.

It seems that there needs to be a splitting and the by using ufl.shape() one can diagnose the issue and find out that v has a 3 dof shape… So I came a little bit further on this and I have a situation that splits and comes through clean without an exception being thrown. I can’t say its perfect yet or without any defect in terms of the math yet but I’ll post this as a tentative solution to the issue of how to assign the forces here in this type of context with 6 dofs and go for a solve:

import dolfinx
from dolfinx import mesh, fem, plot
from mpi4py import MPI
import basix
import numpy as np
from petsc4py import PETSc
import ufl
from dolfinx.fem.petsc import LinearProblem
from dolfinx.fem import dirichletbc, locate_dofs_topological
from ufl import Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, inner, dot
import pyvista as pv

# Material properties
E = 2000000.0  # Young's modulus in PSI (for structural wood)
nu = 0.299999  # Poisson's ratio (for structural wood)
L = 72  # Beam length in inches
W = 4.0  # Beam width in inches
w = 100.0  # Distributed load in lbs/inch

# Define the mesh with tetrahedral cells
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
                         [20, 6, 6], cell_type=mesh.CellType.tetrahedron)

# Define elements for displacements (3 DOF) and rotations (3 DOF)
elem_disp = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(3,))
elem_rot = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(3,))
elem_mixed = basix.ufl.mixed_element([elem_disp, elem_rot])

V = fem.functionspace(domain, elem_mixed)
VV = fem.functionspace(domain, elem_disp)

# Strain tensor as a symmetric gradient of the displacement field
def epsilon(u):
    return 0.5 * (ufl.grad(u) + ufl.grad(u).T)

# Stress tensor derived from Hooke's law
def sigma(u):
    lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))  # First Lamé parameter
    mu = E / (2 * (1 + nu))  # Shear modulus (second Lamé parameter)
    return lambda_ * ufl.tr(epsilon(u)) * ufl.Identity(3) + 2 * mu * epsilon(u)

# Define trial and test functions for displacements (u) and rotations (theta)
(u, theta) = TrialFunctions(V)
(v, tau) = TestFunctions(V)

# Weak form for the elasticity problem (displacement part)
a_disp = inner(sigma(u), epsilon(v)) * ufl.dx

# Rotational-related terms (rotations theta and tau)
a_rot = inner(theta, tau) * ufl.dx

# Coupling terms between displacement and rotation (optional)
a_coupling = inner(theta, v) * ufl.dx + inner(u, tau) * ufl.dx

# Final variational form: combining displacements, rotations, and coupling
a = a_disp + a_rot + a_coupling

# Apply boundary conditions and identify subspaces
V0 = V.sub(0)
V1 = V.sub(1)
Q0, _ = V0.collapse()
Q1, _ = V1.collapse()
fdim = domain.topology.dim - 1

# Locate the entities where x=0 (left boundary) and x=L (right boundary)
def locate_left_boundary(x):
    return np.isclose(x[0], 0.0)

def locate_right_boundary(x):
    return np.isclose(x[0], L)

facets_left = mesh.locate_entities_boundary(domain, fdim, locate_left_boundary)
facets_right = mesh.locate_entities_boundary(domain, fdim, locate_right_boundary)

# Apply boundary condition at the left boundary where x = 0 (clamping)
dofs_left_disp = locate_dofs_topological(V0, fdim, facets_left)
u_left_disp = fem.Function(Q0)
with u_left_disp.vector.localForm() as loc:
    loc.set(0.0)  # Zero displacement in all directions
bc_left_disp = dirichletbc(u_left_disp, dofs_left_disp)

dofs_left_rot = locate_dofs_topological(V1, fdim, facets_left)
u_left_rot = fem.Function(Q1)
with u_left_rot.vector.localForm() as loc:
    loc.set(0.0)  # Zero rotation in all directions
bc_left_rot = dirichletbc(u_left_rot, dofs_left_rot)

# Apply boundary condition at the right boundary where x = L (clamping)
dofs_right_disp = locate_dofs_topological(V0, fdim, facets_right)
u_right_disp = fem.Function(Q0)
with u_right_disp.vector.localForm() as loc:
    loc.set(0.0)  # Zero displacement in all directions
bc_right_disp = dirichletbc(u_right_disp, dofs_right_disp)

dofs_right_rot = locate_dofs_topological(V1, fdim, facets_right)
u_right_rot = fem.Function(Q1)
with u_right_rot.vector.localForm() as loc:
    loc.set(0.0)  # Zero rotation in all directions
bc_right_rot = dirichletbc(u_right_rot, dofs_right_rot)

# Combine all boundary conditions
bcs = [bc_left_disp, bc_left_rot, bc_right_disp, bc_right_rot]

# Create functions for displacements and rotations
f_disp = fem.Function(Q0)
f_rot = fem.Function(Q1)

# Compute the distributed load per node and apply it
nodes = domain.geometry.x
num_nodes = len(nodes)

# Assign distributed force (downward in -y) and moment (around x-axis)
I = (W**4) / 12  # Moment of inertia for a rectangular cross-section

# Loop over all nodes to apply the distributed load and moment
for i, node in enumerate(nodes):
    fy = -w * L / num_nodes  # Distributed load down the -y axis
    y = node[1]  # y-coordinate of the node
    mx = -(w * abs(y) * L**2) / (2 * I)  # Moment around x-axis

    # Apply fy to the displacement component in y (index 1)
    f_disp.vector.setValue(i * 3 + 1, fy)  # Displacement force in y

    # Apply mx to the rotational component around x (index 0)
    f_rot.vector.setValue(i * 3 + 0, mx)  # Rotational moment around x

# Combine the force vectors manually
f = fem.Function(V)
with f.vector.localForm() as f_local, f_disp.vector.localForm() as disp_local, f_rot.vector.localForm() as rot_local:
    f_local.set(0)  # Initialize the combined force vector to zero
    indices_disp = np.arange(num_nodes * 3, dtype=np.int32)  # Indices for displacement components
    indices_rot = np.arange(num_nodes * 3, dtype=np.int32)  # Indices for rotational components
    f_local.setValuesLocal(indices_disp, disp_local[:])
    f_local.setValuesLocal(indices_rot + 3, rot_local[:])  # Offset for rotational components
    f_local.assemble()  # Finalize the combined force vector


# Define the right-hand side (RHS) for the variational problem
l = inner(f.sub(0), v) * ufl.dx + inner(f.sub(1), tau) * ufl.dx  # Adjust RHS to match the shape of the test functions

# Now create the linear problem
problem = LinearProblem(a, l, bcs=bcs)
uh = problem.solve()

# Save the solution to a CSV file
np.savetxt('U.csv', uh.x.array)
print('done...')

array([-4.53852396e-03,  3.19832883e-04,  5.44793081e-04, ...,
       -7.86804411e-05, -4.22010879e-07, -4.28360461e-07])

If there is any comments further I would very much like to hear them, that is why I won’t be marking a solution at this time.

So at this point I can clearly see that there is a problem in accurate depiction so far:

So I don’t know at this time what went wrong, however I would suspect at this point that the force vector’s locations have a winding issue to deal with… For sure this looks a no good depiction of the story so far that should be coming up…

Any help about how to clear that up would be greatly appreciated!

[prusso630]

Oh… I guess this is what that is being said here:
https://bleyerj.github.io/comet-fenicsx/intro/plates/plates.html

It is just now coming to my attention, so the intro says the classical hermits mockup isn’t really covered here so there is an alternate way to look at in the docs… So its more advanced, so far it looks like the moments needed are applied with gmsh…

Although the code that was posted is a way to start a mixed element scenario it looks like the ufl is no good there… not at sure if such a way of doing things can ever work maybe that is what is being covered in the 3d beam example shown in the link…

[prusso630]

In my efforts to do this I got really confused and into some disarray… I was a bit absent minded it seems that this type of problem can maybe be constructed as a hermite system with fenics but not yet with dolfinx… So once that had dawned on me after getting this mixed element problem going, I switched over to another way of doing this:

How to solve the 3D beam problem for a rectangular beam of linear tetrahedral elements

dolfinx version: 0.8.0 Here is a script to generate the beam: import numpy as np import gmsh import ufl from mpi4py import MPI from dolfinx import mesh, fem, io import basix from petsc4py import PETSc import dolfinx.fem.petsc # Initialize gmsh and load your model gmsh.initialize() gmsh.open("beam_mesh.msh") # Create the DOLFIN-X mesh using gmshio.model_to_mesh domain, _, _ = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0) # Finalize gmsh gmsh.finalize() gdim = domain.geometry.dim t…

As far as why are the results here no good… somewhere in discussion is a way to hermits with fenics that looks promising if one wants it this way or similar…