Trying to have cantilever beam case Converting legacy FEniCS code to FEniCSx to work, I want to plot displacement to make sure loads and boundaries are ok.
I fetch DOFs coordinates from :
dof_coords = W0.tabulate_dof_coordinates()
… and y displacements from :
dof_displacements_y = u.sub(0).sub(1).collapse().x.array
which gives me this sawtooth curve and raises interrogations on index discrepancy between the values :
I have been through different threads including Recover deformation field vectors from solution but could not find how to fix my case.
I would greatly appreciate any insights or suggestions on how to resolve this issue.
Here bellow my full code :
import os
arch = os.getenv("ARGS", "real")
try:
import google.colab # noqa: F401
except ImportError:
import ufl
import dolfinx
else:
try:
import ufl
import dolfinx
except ImportError:
if arch != "complex":
!wget "https://fem-on-colab.github.io/releases/fenicsx-install-real.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh"
else:
!wget "https://fem-on-colab.github.io/releases/fenicsx-install-complex.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh"
import ufl
import dolfinx
# Install meshio for mesh handling
!pip install meshio
# Importing necessary libraries and modules
import numpy as np # Numerical operations
from mpi4py import MPI # Parallel computing
import basix # Finite element basis functions
from dolfinx import fem
from dolfinx import plot
from dolfinx.fem import (Constant, FunctionSpace, locate_dofs_geometrical,
Function, dirichletbc, assemble_scalar, form)
# FEM module functions
from dolfinx.io import XDMFFile # I/O operations
from dolfinx.fem.petsc import LinearProblem
# Solving linear problems
from dolfinx.mesh import locate_entities, meshtags
# Mesh utility functions
from ufl import (Jacobian, as_vector, sqrt, inner, cross, dot, TestFunction,
TrialFunction, split, diag, grad, Measure, dx, avg)
# Defining variational forms
import meshio # Mesh I/O operations
import matplotlib.pyplot as plt
import pyvista as pv # Visualization
# Initialize PyVista
pv.start_xvfb()
## Simple mesh
length = 10.0
N = 10 # number of elements
points = np.zeros((N+1, 3)) # N+1 array of zeros in 3D space
points[:, 0] = np.linspace(0, length, N+1) # Distribute points along x from 0 to 'length'
cells = np.array([[i, i+1] for i in range(N)]) # line elements connectivity
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells})) # Write XDMF mesh file
# Reading the mesh from the XDMF file
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid")
# Creating connectivity information for the mesh
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)
# Explicitly compute the entity-to-cell connectivity
# domain.topology.create_connectivity(domain.topology.dim, 0)
# mesh test
# Get the geometrical and topological dimensions of the mesh from the 'domain' object.
gdim = domain.geometry.dim # = 3
tdim = domain.topology.dim # = 1
# Print the geometrical and topological dimensions of the mesh.
print(f"Geometrical dimension = {gdim}")
print(f"Topological dimension = {tdim}")
# Extract coordinates of vertices
vertices_coordinates = domain.geometry.x
print("len(vertices_coordinates)=",len(vertices_coordinates)) # = 11
# Create a 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Plot vertices
ax.scatter(vertices_coordinates[:, 0],
vertices_coordinates[:, 1],
vertices_coordinates[:, 2],
c='b',
marker='o')
# Set labels
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
# Show plot
plt.show()
# Defining a function to compute the tangent vector to the mesh
def tangent(domain):
t = Jacobian(domain)
return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))
# Calculating the tangent vector for the domain
t = tangent(domain)
# Defining basis vectors orthogonal to the tangent vector
ez = as_vector([0, 0, 1]) # Basis vector in the z-direction
a1 = cross(t, ez) # First orthogonal vector
a1 /= sqrt(dot(a1, a1)) # Normalizing the first orthogonal vector
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))
# Geometric properties
thick_val = 0.3
thick = Constant(domain,thick_val)
width_val = thick_val
width = thick/3
# Material Properties
E_val = 70e3
E = Constant(domain,E_val) # Young's modulus
nu = Constant(domain,0.3) # Poisson's ratio
G = E/2/(1+nu) # Shear modulus
rho = Constant(domain,2.7e-3) # Density
g = Constant(domain,9.81) # Gravitational acceleration
S = thick*width # Cross-sectional area
ES = E*S # Axial stiffness
EI1 = E*width*thick**3/12 # Bending stiffness about one axis
EI2 = E*width**3*thick/12 # Bending stiffness about the other axis
J = (width * thick**3 + thick * width**3) / 3
GJ = G*J # Torsional stiffness
kappa = Constant(domain,5./6.) # Shear correction factor
GS1 = kappa*G*S # Shear stiffness in one direction
GS2 = kappa*G*S # Shear stiffness in the other direction
# Define the displacement vector element using basix
Ue = basix.ufl.element(
family=basix.ElementFamily.P, # Lagrange family
cell=basix.CellType.interval, # 1D cell type
degree=2, # Polynomial degree
shape=(3,)) # 3 components vector
# Define the rotation vector element using basix
Te = basix.ufl.element(
family=basix.ElementFamily.P, # Lagrange family
cell=basix.CellType.interval, # 1D cell type
degree=1, # Polynomial degree
shape=(3,)) # 3 components vector
# Creating a mixed element for displacement and rotation
Ueprod = basix.ufl.mixed_element([Ue, Te]) # Mixed element
W = dolfinx.fem.functionspace(domain, Ueprod) # Creating a function space
# Defining test and trial functions
u_ = TestFunction(W) # Test function
du = TrialFunction(W) # Trial function
(w_, theta_) = split(u_) # Split test function into displacement and rotation
(dw, dtheta) = split(du) # Split trial function into displacement and rotation
# Defining generalized strains function
def tgrad(u):
return dot(grad(u), t)
def generalized_strains(u):
(w, theta) = split(u)
return as_vector([dot(tgrad(w), t),
dot(tgrad(w), a1)-dot(theta, a2),
dot(tgrad(w), a2)+dot(theta, a1),
dot(tgrad(theta), t),
dot(tgrad(theta), a1),
dot(tgrad(theta), a2)])
# Defining generalized stresses function
def generalized_stresses(u):
return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])),
generalized_strains(u))
# Calculating stresses and strains
Sig = generalized_stresses(du) # Generalized stresses
Eps = generalized_strains(u_) # Generalized strains
# Locate left boundary
def left(x):
return(np.isclose(x[0], 0, 0.1)) # value, target, tolerance
# Locate load application
def load(x):
return np.isclose(x[0], length, 0.1)
# Defining boundaries and marking facets
boundaries = [(1, left), # Marker 1 for left boundary
(2, load)] # Marker 2 for load application point
facet_indices, facet_markers = [], [] # Lists to store facet indices and markers
fdim = domain.topology.dim - 1 # Facet dimension
# Locating and marking facets for each boundary
for (marker, locator) in boundaries:
facets = locate_entities(domain, fdim, locator) # Locate facets
facet_indices.append(facets) # Add facet indices
facet_markers.append(np.full_like(facets, marker)) # Add facet markers
facet_indices = np.hstack(facet_indices).astype(np.int32)# Combine facet indices
facet_markers = np.hstack(facet_markers).astype(np.int32)# Combine facet markers
sorted_facets = np.argsort(facet_indices) # Sort facets
facet_tag = meshtags(domain,
fdim,
facet_indices[sorted_facets],
facet_markers[sorted_facets]) # Create mesh tags
# Collapsing subspaces of the function space
W0, submap0 = W.sub(0).collapse() # Collapsing displacement field subspace
W1, submap1 = W.sub(1).collapse() # Collapsing rotation field subspace
left_dofs_x = locate_dofs_geometrical((W.sub(0),W0),left )
left_dofs_r = locate_dofs_geometrical((W.sub(1),W1),left)
# Create zero function for the first subspace (W0)
zero_0 = Function(W0)
with zero_0.vector.localForm() as zero_local:
zero_local.set(0.0)
# Create zero function for the second subspace (W1)
zero_1 = Function(W1)
with zero_1.vector.localForm() as zero_local:
zero_local.set(0.0)
# Define boundary conditions
bcs = [
dirichletbc(zero_0, left_dofs_x, W.sub(0)),
dirichletbc(zero_1, left_dofs_r, W.sub(1))]
metadata = {"quadrature_scheme": "default", "quadrature_degree": 4}
# Define a measure for surface integrals
ds = Measure("ds",
domain=domain,
subdomain_data=facet_tag,
metadata = {"quadrature_degree": 4})
# Define a measure for surface integrals in 3D
dS = Measure("dS",
domain=domain,
subdomain_data=facet_tag,
metadata = {"quadrature_degree": 4})
# Define a measure for volume integrals, used for shear stress
dx_shear = dx(scheme="default",
metadata=metadata)
# Define the bilinear form (stiffness matrix)
k_form = (sum([Sig[i]*Eps[i]*dx for i in [0, 3, 4, 5]]) # normal and bending terms
+(Sig[1]*Eps[1]+Sig[2]*Eps[2]) # shear terms
*dx_shear) # Integrate
P = -1.0 # load intensity
# Define the linear form (load vector)
l_form = (inner(Constant(domain,P), # apply P
avg(u_[1])) # along y
*dS(2)) # on mark 2
u = Function(W) # Define the solution function
# Define linear problem
problem = LinearProblem(k_form,
l_form,
u=u,
bcs=bcs)
uh=problem.solve() # Solve the linear problem
xyz = W0.tabulate_dof_coordinates()
x_coor = xyz[:,0]
y_coor = xyz[:,1]
z_coor = xyz[:,2]
x_value, y_value, z_value = [], [], []
for i in range(len(x_coor)):
x_value.append(u.sub(0).x.array[3*i:3*(i+1)][0])
y_value.append(u.sub(0).x.array[3*i:3*(i+1)][1])
z_value.append(u.sub(0).x.array[3*i:3*(i+1)][2])
# Get the coordinates of the DoFs in the collapsed subspace
dof_coords = W0.tabulate_dof_coordinates()
# Get the y-displacements in the collapsed subspace
dof_displacements_y = u.sub(0).sub(1).collapse().x.array
# Ensure `submap0` does not contain out-of-bound indices
valid_indices = [i for i in submap0 if i < len(dof_displacements_y)]
# Initialize mapped_displacements array with zeros
mapped_displacements = np.zeros(len(valid_indices))
# Map the displacements
for i, original_index in enumerate(valid_indices):
mapped_displacements[i] = dof_displacements_y[original_index]
# Now you can plot or use `mapped_displacements` with `dof_coords`
plt.plot(dof_coords[:, 0], mapped_displacements, 'o-')
plt.xlabel('X-coordinate')
plt.ylabel('Y-displacement')
plt.title('Y-displacements along the beam')
plt.show()
# Print the maximum y-displacement
print("Max y-displacement =", min(mapped_displacements))
# Print the maximum y-displacement
dof_displacements_y = u.sub(0).sub(1).collapse().x.array
print("dof_displacements_y = ",dof_displacements_y)
print ("len(dof_displacements_y) = " ,len(dof_displacements_y)) # = 21 because degree=2
# Get DOFs coordinates for the displacement field
dof_coords = W.sub(0).collapse()[0].tabulate_dof_coordinates()
# Plot the data
plt.plot(dof_coords[:, 0], dof_displacements_y)
# Add labels and title
plt.xlabel('x_coordinates')
plt.ylabel('y_displacements')
plt.title('DOFs')
# Show the plot
plt.show()
print("Max y-displacement =",
min(dof_displacements_y)) # = -161.5063564124886
# Calculate maximum y-displacement using beam theory
beam_theory_deflection = (P * length**3) / (3 * E_val*width_val*thick_val**3/12)
print("Beam theory =", beam_theory_deflection) # = -7.054673721340389
