I modified the code a bit. Now it does run and the spring is correctly setup at the first and last node but not in nodes in the middle. What am I doing wrong?
Part of the code I modified:
domain.topology.create_connectivity(0, domain.topology.dim)
pts = np.array([], dtype=np.int32)
for pt in supports:
startpt=locate_entities(domain,0,lambda x : np.isclose(x[0], pt))
pts = np.concatenate((pts, startpt))
tag = meshtags(domain, int(domain.topology.dim - 1), pts, np.ones_like(pts))
ds = Measure("ds", domain=domain, subdomain_data=tag)
k_form += k_spring * inner(v, u_) * ds(1)
Full code:
import numpy as np
from mpi4py import MPI
import basix
from dolfinx.fem import (Function, dirichletbc,functionspace, form, set_bc, assemble_scalar,
locate_dofs_topological,Constant,Expression)
from dolfinx.io import VTKFile
from dolfinx.fem.petsc import LinearProblem, assemble_matrix, assemble_vector, apply_lifting
from dolfinx.mesh import locate_entities, locate_entities_boundary,create_interval, meshtags
from ufl import (SpatialCoordinate,inner, dot, TestFunction, TrialFunction, action, div, grad, dx, Measure, derivative, conditional, le, lhs, rhs)
import matplotlib.pyplot as plt
from section import inertia, area
import dolfinx.fem.petsc as petsc
import time
t1 = time.time()
def last_value(arr):
if len(arr) == 0:
return -1
return arr[-1]
#################################################################
########### SUPPORTS PARAMETERS #################################
#################################################################
supports = [0.0,5.0,10.0]
k_spring = 4000 / 0.004
########## GEOMETRIC INPUT ####################
E = 210e9 # N/m^2
L = 10.0 # m
thickness = 0.005 # m
h = 0.12 # m
b = 0.09 # m
rho= 7.8e-3
g = 9.81
#################################################################
########### CONSTRUCT BEAM MESH #################################
#################################################################
NUM_ELEM = int(L * 100)
domain = create_interval(MPI.COMM_WORLD, NUM_ELEM, [0, L])
#################################################################
##### ENTER MATERIAL PARAMETERS AND CONSTITUTIVE MODEL ##########
#################################################################
x = SpatialCoordinate(domain)
E = Constant(domain,E)
rho = Constant(domain,rho)
g = Constant(domain,g)
k_spring = Constant(domain,k_spring)
f = Constant(domain,1000.0)
I = b * h ** 3 / 12
EI = E * I
#################################################################
########### COMPUTE STATIC SOLUTION #############################
#################################################################
#define Moment expression
def M(u):
return EI*div(grad(u))
# Create Hermite order 3 on a interval (for more informations see:
# https://defelement.com/elements/examples/interval-Hermite-3.html )
beam_element=basix.ufl.element(basix.ElementFamily.Hermite, basix.CellType.interval, 3)
#finite element function space on domain, with trial and test fxns
W = functionspace(domain,beam_element)
print("Number of DOFs: %d" % W.dofmap.index_map.size_global)
print("Number of elements (intervals): %d" % NUM_ELEM)
print("Number of nodes: %d" % (NUM_ELEM+1))
u_ = TestFunction(W)
v = TrialFunction(W)
#bilinear form (LHS)
k_form = inner(div(grad(v)),M(u_))*dx
#linear form construction (RHS)
l_form = f*u_*dx
#APPLY BOUNDARY CONDITIONS
#initialize function for boundary condition application
ubc = Function(W)
with ubc.vector.localForm() as uloc:
uloc.set(0.)
# list of boundary conditions
bcs = []
domain.topology.create_connectivity(0, domain.topology.dim)
pts = np.array([], dtype=np.int32)
for pt in supports:
startpt=locate_entities(domain,0,lambda x : np.isclose(x[0], pt))
pts = np.concatenate((pts, startpt))
tag = meshtags(domain, int(domain.topology.dim - 1), pts, np.ones_like(pts))
ds = Measure("ds", domain=domain, subdomain_data=tag)
k_form += k_spring * inner(v, u_) * ds(1)
#SOLVE VARIATIONAL PROBLEM
#initialize function in functionspace for beam properties
u = Function(W)
# solve variational problem
problem = LinearProblem(k_form, l_form, u=u, bcs=bcs)
uh=problem.solve()
#uh.name = "Displacement and Rotation"
uh_grad = grad(uh)
uh_xx = grad(uh_grad[0])
uh_xxx = grad(uh_xx[0])
# Define the shear force expression
shear_force_expr = EI * uh_xxx
V = functionspace(domain, ("DG", 0))
shear_expr = Expression(shear_force_expr, V.element.interpolation_points())
shear = Function(V)
shear.interpolate(shear_expr)
reaction = []
for support in supports:
index=locate_entities(domain,0,lambda x : np.isclose(x[0], support))[0]
if index == 0:
left_side = 0
else:
left_side = shear.x.array[index-1]
if index == len(shear.x.array):
right_side = 0
else:
right_side = shear.x.array[index]
f_reac = right_side - left_side
reaction.append(f_reac)
print(f"Support at {support} has reaction force of {round(f_reac / 1000, 2)} kN")
print("time taken: ", time.time()-t1)
#################################################################
########### SAVE AND VISUALIZE RESULTS ##########################
#################################################################
#NOTE: The solution uh contains both the rotation and the displacement solutions
#The rotation and displacment solutions can be separated as follows:
disp = np.empty(0)
rot = np.empty(0)
for i,x in enumerate(uh.x.array):
if i % 2 != 0:
rot = np.append(rot,x)
else:
disp = np.append(disp,x)
x_fem = np.linspace(0,L,disp.shape[0])
figure, axis = plt.subplots(1,1)
print("Maximum magnitude displacement is: %e" % np.min(disp))