How to obtain a solution for a problem with 6 degrees of freedom per element?

[prusso630]

print(f'dolfinx: {dolfinx.__version__}')
'''
dolfinx: 0.8.0
'''

import dolfinx
from dolfinx import mesh, fem
from mpi4py import MPI
import basix
import numpy as np
from petsc4py import PETSc

# A six degrees of freedom per node tetrahedral domain was created and a #function space was created with the element called V...
...


# Create a function to hold the force vector at each node
force = fem.Function(V)

# Initialize the force vector to zero
force.vector.set(0.0)

# Get the topology and connectivity of the mesh
topology = domain.topology 
cells = topology.connectivity(topology.dim, 0)
cell_vertices = domain.geometry.x[cells.array]

# Define a function to calculate the force vector at each node based on the coordinates
def calculate_force(x):
    fx = 0.1 * x[0]  # Example formula for force in x-direction
    fy = 0.2 * x[1]  # Example formula for force in y-direction
    fz = 0.3 * x[2]  # Example formula for force in z-direction
    mx = 0.0         # Example moment about x-axis (set to 0)
    my = 0.0         # Example moment about y-axis (set to 0)
    mz = 0.0         # Example moment about z-axis (set to 0)
    return [fx, fy, fz, mx, my, mz]

# Assemble the force vector by looping over the DOFs
dofmap = V.dofmap
num_dofs = force.vector.getSize()

# Define the bilinear form (stiffness matrix)
u = fem.Function(V)
v = fem.Function(V)

# Access the vector of the function
u_vector = u.vector.getArray()
v_vector = u.vector.getArray()
# Extract displacements (assuming the first 3 DOFs correspond to displacements)
u_disp = u_vector[:3]  # Get the first 3 components for displacement (u1, u2, u3)
v_disp = v_vector[:3]
# Create UFL expressions for the displacement components (use vector components directly)
u_disp = u.split()[0]  # Extract the displacement component (u1, u2, u3)
v_disp = v.split()[0]  # Extract the displacement component (v1, v2, v3)

# Compute the gradient of each displacement component
grad_u_disp = ufl.grad(u_disp)  # Gradient of displacement u (which will give us a 3x3 tensor)
grad_v_disp = ufl.grad(v_disp)  # Gradient of displacement v (same 3x3 tensor)

# Now compute the symmetric part of the gradient (this is the strain tensor)
epsilon_u = ufl.sym(grad_u_disp)  # Symmetric part of the gradient for strain tensor of u
epsilon_v = ufl.sym(grad_v_disp)  # Symmetric part of the gradient for strain tensor of v



# For rotations, define an appropriate formulation (e.g., assuming small rotations)
# You may need to modify this based on your problem's specific formulation.
# For example, rotation-related strain could be modeled in a different way.
# Here, we assume the rotation strain is zero for simplicity.
#epsilon_rot_u = u[3:6]  # Extract rotations (last 3 DOFs)
#epsilon_rot_v = v[3:6]  # Extract rotations (last 3 DOFs)

# Bilinear form for displacements
#a_displacements = fem.form(ufl.inner(epsilon_u, epsilon_v) * ufl.dx)

# Bilinear form for rotations (example, can be expanded for real problem)
#a_rotations = fem.form(ufl.inner(epsilon_rot_u, epsilon_rot_v) * ufl.dx)

# Total bilinear form (sum of displacement and rotation contributions)
#a = a_displacements + a_rotations

# Define the linear form (force vector)
#L = fem.form(force * v)  # Linear term (force)

# Create a linear problem
#problem = fem.LinearProblem(a, L, petsc_options={"ksp_type": "gmres", "pc_type": "hypre"})

# Solve the linear problem
#solution = problem.solve()

# Optionally, print out the solution (displacement and rotations) at nodes
#print(solution.vector)

So far I am under the impression that a more advanced type of bi-linear form would have to be created to solve a problem with forces and rotations. This is my first attempt to do so. I’ll admit that so far I am just understanding what the trial and test function actually are and I am working under the impression that u and v stem from gradients that are formed from the evaluation of the basis functions and some other factors during assembly. I have a text that describes the process fairly well and some MatLabs that I read through…

At this point, I can see that epsilon_u = ufl.sym(grad_u_disp), is certainly not the way to be doing things and at this point does lead to:

Traceback (most recent call last):
  File "/home/prusso/dolfinx-beam_w/tetra_forcevector3.py", line 66, in <module>
    epsilon_u = ufl.sym(grad_u_disp)  # Symmetric part of the gradient for strain tensor of u
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/operators.py", line 337, in sym
    return Sym(A)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/tensoralgebra.py", line 542, in __new__
    raise ValueError("Symmetric part of tensor with rank != 2 is undefined.")
ValueError: Symmetric part of tensor with rank != 2 is undefined.

So if I understand correctly the element here of degree 1 can be solved using “LinearProblem”…

Could I get some help to get a solution vector U coming out here? One thing I am also inquiring is, how to set a demonstration force vector. Mostly I would like to say for example have forces on the Y-axis only of the same amount each node. I will have to set up Dirichlet conditions as well and I understand that those aren’t shown since this is a work in progress still…I plan to have a fixed support so far at all nodes x=0.0…

Is anything on the right track here? How to start to get a solution vector U from this type of problem?

[prusso630]

As an aside… I made an attempt to set up a Dirichlet condition for this 6 degrees of freedom per node problem and I ran into a shape and size mismatch that I am not really sure at this point what to do with:

# Define the boundary condition value (zero displacement and zero rotation)
bc_value = fem.Function(V)
bc_value.interpolate(lambda x: np.zeros((6,)))

# Locate nodes where x = 0
def boundary_x0(x):
    return np.isclose(x[0], 0.0)

# Apply Dirichlet boundary conditions to these nodes
bc = fem.dirichletbc(bc_value, fem.locate_dofs_geometrical(V, boundary_x0), V)

the difficulty seems to occur at the interpolation lambda… I guess I will also have to work through this otherwise I won’t be able to determine a solution vector…

So far as far as setting up a force vector this is the best I have where element is a 6 d.o.f. per node tetrahedron in this case…

# Create the function space with the defined element
V = fem.functionspace(domain, element)

# Create a function to hold the force vector at each node
force = fem.Function(V)

# Initialize the force vector
force.vector.set(0.3)

I guess this is maybe the cause to split the forces and rotations into two different elements, so in my efforts to come to a solution I was thrown off track a bit… So far I am a bit confused still how to come to a congruent solution from all 6 dofs… It seems there is more to this than I really fully understand so far… If the idea is to use two separate 3 dof elements could I see som code to that effect ho to come to a congruent solution… So far the rotational forces influence the transvers forces and vice versa for all I know…

So if the translational and rotational forces are compartmentalized things need to come back together sometime at solution time like a swing back together of the bilinear form. I have an idea how that might work however I prefer to ask the question to prevent the making of an assumption here since that it would be very prohibitive to have to test all of these things myself…

So far from what I know commonly there is B matrix of partial derivatives that leads to the influence of all degrees of freedom. So in this case for LinearSolver there is the input of some combination of the bilinear form and the test or trial function and the petsc maybe uses the right hand vector although that is a blackbox to me so far…

So if it is known how to do this I would really appreciate a reply…

Also, if someone can point me in the direction if there is code to look at how the Linear solver makes use of the input parameters that would also be great!

[prusso630]

I got a little bit further with this in terms of getting to the point where some type of solving of the problem can maybe go on… Not sure if this will lead to a solution that is good yet. Here is the code I have:

import dolfinx
from dolfinx import mesh, fem
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_geometrical


L = 72
W = 4.0

# 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 and rotations
element_disp = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(3,))
element_rot = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(3,))

# Create function spaces for displacements and rotations
V_disp = fem.functionspace(domain, element_disp)
V_rot = fem.functionspace(domain, element_rot)

# Create functions to hold the force vectors at each node
force_disp = fem.Function(V_disp)
force_rot = fem.Function(V_rot)

# Initialize the force vectors to zero
force_disp.vector.set(0.3)
force_rot.vector.set(0.0)

# Define functions for displacement and rotation
u_disp = fem.Function(V_disp)
u_rot = fem.Function(V_rot)
v_disp = fem.Function(V_disp)
v_rot = fem.Function(V_rot)

# Define UFL expressions for the displacement and rotation components
grad_u_disp = ufl.grad(u_disp)  # Gradient of displacement u
grad_v_disp = ufl.grad(v_disp)  # Gradient of displacement v

# Symmetric part of the gradient (strain tensor) for displacements
epsilon_u = ufl.sym(grad_u_disp)
epsilon_v = ufl.sym(grad_v_disp)

# Bilinear form for displacements
a_disp = fem.form(ufl.inner(epsilon_u, epsilon_v) * ufl.dx)

# (Example) Bilinear form for rotations, assuming small rotations
a_rot = fem.form(ufl.inner(u_rot, v_rot) * ufl.dx)

# Total bilinear form (sum of displacement and rotation contributions)
a = fem.form(ufl.inner(epsilon_u, epsilon_v) * ufl.dx + ufl.inner(u_rot, v_rot) * ufl.dx)



# Define the linear form (force vector)
L_disp = fem.form(ufl.dot(force_disp, v_disp) * ufl.dx)
L_rot = fem.form(ufl.dot(force_rot, v_rot) * ufl.dx)

# Define the total linear form
L = fem.form(ufl.dot(force_disp, v_disp) * ufl.dx + ufl.dot(force_rot, v_rot) * ufl.dx)



# Define the boundary condition
def boundary(x):
    return np.isclose(x[0], 0)

# Define the Dirichlet boundary value for vector-valued spaces
zero_disp = fem.Constant(domain, PETSc.ScalarType((0.0, 0.0, 0.0)))
zero_rot = fem.Constant(domain, PETSc.ScalarType((0.0, 0.0, 0.0)))

# Apply Dirichlet boundary conditions to the displacements and rotations
bc_disp = dirichletbc(zero_disp, locate_dofs_geometrical(V_disp, boundary), V_disp)
bc_rot = dirichletbc(zero_rot, locate_dofs_geometrical(V_rot, boundary), V_rot)




# Combine boundary conditions
bcs = [bc_disp, bc_rot]

# Define the linear problem
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

# Solve the linear problem
u_sol = problem.solve()

# Split the solution into displacement and rotation components
u_disp, u_rot = u_sol.split()

# Output the results
print(f"Displacement:\n{u_disp.x.array}")
print(f"Rotation:\n{u_rot.x.array}")

print(f'dolfinx version: {dolfinx.__version__}')

The difficulty I have is:

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'

So there seems to be a problem with using LinearProblem this way. Can I get some more information on what is causing this breakdown in the solver and if I am on the right track to obtaining a solution here? The reason for my edit was the lack of Dirichlet here which are now included some way…

[prusso630]

I finally reached the understanding that using a 6dof element was the cause of the dilemma. Later on it was found out how to handle this type of formulation by splitting it into a mixed element situation:

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

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…

and also through all this I became a bit absent minded. So it was brought to my attention some time ago that the classical hermits formulation of six dofs can’t work like this with dolfinx but there could be a way to do so with fenics and also there could be support later on based from what I read for a hermife style of formulation.