Targeting subspaces within MixedElement function spaces

1

Hi,

I’m currently writing a script to calculate the deflection and twist of a 1D beam. This beam is horizontal up until a point where it has a dihedral section at a 25 degree angle. I have been able to apply a vertical force to this using PointSource and targeting the subspace of my MixedElement function space. My question is, how would I go about applying a force in the horizontal direction? I’ve done this before in 2D without a MixedElement function space by targeting the subspace of a vector function space but I’m unsure how to do this with a MixedElement function space.

Here’s my code:

import meshio
from matplotlib import pyplot as plt
import time
from dolfin import *

t0 = time.time()

# Read in mesh
msh = meshio.read("Linemesh.msh")
meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"line": msh.cells["line"]}))
mesh = Mesh()
filename = "mesh.xdmf"
f = XDMFFile(filename)
f.read(mesh)
k = 1

# Define function space
UE = FiniteElement("CG",mesh.ufl_cell(),k) # Displacement u
VE = FiniteElement("CG",mesh.ufl_cell(),k) # EI*u''
Phi = FunctionSpace(mesh,"CG", k) # Phi
W = FunctionSpace(mesh,MixedElement([UE,VE]))

# This is the most versatile way to use spatial coordinates in a formulation:
x = SpatialCoordinate(mesh)

# Characteristic functions for each end of the beam, for use in boundary
# integrals for weak enforcement of boundary conditions:
leftChar = 10000-x[0]
rightChar = x[0]

# Defining properties
# Outer diameters
D0 = 30 #mm
D1 = 30 #mm
D2 = 30 #mm
Do = Expression('x[0] <= 250 ? D0 : ((250 < x[0] && (x[0] <= 500)) ? D1 : D2)'\
               ,degree = 0, D0 = D0, D1 = D1, D2 = D2) #mm

# Inner diameters
D0 = 25 #mm
D1 = 25 #mm
D2 = 25 #mm
Di = Expression('x[0] <= 250 ? D0 : ((250 < x[0] && (x[0] <= 500)) ? D1 : D2)'\
               ,degree = 0, D0 = D0, D1 = D1, D2 = D2) #mm

I = (3.141592654/64)*((Do**4)-(Di**4)) #mm^4

# Young's modulus
E0 = 150*10**3 
E1 = 150*10**3 
E2 = 150*10**3 
E = Expression('x[0] <= 250 ? E0 : ((250 < x[0] && (x[0] <= 500)) ? E1 : E2)'\
               ,degree = 0, E0 = E0, E1 = E1, E2 = E2) #Pa   
EI = E*I
J = (3.141592654/32)*((Do**4)-(Di**4))#mm^4

# Shear modulus
G0 = 75*10**3
G1 = 75*10**3 
G2 = 75*10**3 
G = Expression('x[0] <= 250 ? G0 : ((250 < x[0] && (x[0] <= 500)) ? G1 : G2)'\
               ,degree = 0, G0 = G0, G1 = G1, G2 = G2) #Pa   
GJ = G*J
qx = Constant(0) #Nm/m Uniformly distributed torque
q = Constant(0) #N/m Uniformly distributed lift

n = FacetNormal(mesh)
h = CellDiameter(mesh)

# PDE residual; BC enforcement not immediately obvious!
def pdeRes(u,v,du,dv):
    return inner(grad(u),grad(du))*dx + inner(v/EI,du)*dx \
        + inner(grad(v),grad(dv))*dx + inner(q,dv)*dx \
        + leftChar*u*inner(n,grad(dv))*ds \
        + rightChar*v*inner(n,grad(du))*ds \
        - leftChar*dot(grad(v),n)*dv*ds \
        - rightChar*dot(grad(u),n)*du*ds

# Set up functions and assemble
u,v= split(TrialFunction(W))
du,dv = split(TestFunction(W))
residual = pdeRes(u,v,du,dv)
a, L = lhs(residual), rhs(residual)
w = Function(W)
A = assemble(a)
b = assemble(L)

# Torque forumulation
phi = TrialFunction(Phi)
dphi = TestFunction(Phi)
LT = qx*dphi*dx # Linear part
aT = inner(GJ*grad(phi),grad(dphi))*dx #Bilinear part

# Set BCs and solve:
bc = DirichletBC(Phi,Constant(0.0),"x[0]<DOLFIN_EPS")
AT = assemble(aT)
bT = assemble(LT)

# Apply PointSources for lift
x = 10000
y = 0
z = 1056.546
sub1 = W.sub(1)
ptSrc = PointSource(sub1, Point(x, y ,z), 10)
ptSrc.apply(b)

# Apply PointSources for conentrated torque (Nmm)
ptSrc1 = PointSource(Phi, Point(6660,0), 0.01)
ptSrc2 = PointSource(Phi, Point(0,0), 0)
ptSrc3 = PointSource(Phi, Point(0, 0), 0)
ptSrc1.apply(bT)
ptSrc2.apply(bT)
ptSrc3.apply(bT)

# Solve
u,v = w.split()
phi_h = Function(Phi)
solve(A, u.vector() ,b)
solve(AT ,phi_h.vector(),bT)

# Plot deflection:
#plot(u)
#ax1 = plt.subplots(1)
                    
V = FunctionSpace(mesh, 'CG', 2)
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
print('max deflection:',
u_magnitude.vector().get_local().max())

# Compute magnitude of displacement
V = FunctionSpace(mesh, 'CG', 2)
u_magnitude = sqrt(dot(phi_h, phi_h))
u_magnitude = project(u_magnitude, V)
print('max twist:',
u_magnitude.vector().get_local().max())

fileResults = XDMFFile('Results.xdmf')
fileResults.parameters["flush_output"] = True
fileResults.parameters["functions_share_mesh"] = True

U = FunctionSpace(mesh, UE)
u = project(u, U)
disp = Function(U, name='Displacement')
disp.assign(u)
fileResults.write(disp,0.)

Twist = Function(V, name='Twist')
Twist.assign(phi_h)
fileResults.write(Twist,0.)

fileResults.close()
t1 = time.time()

total = t1-t0
print(total)

I’m guessing that using FiniteElement for UE might be incorrect if I need to apply forces at two different vectors, but I’m unsure what the alternative would be.

Any help would be much appreciated. Thanks.

Edit: Sorry about the long code, the important bits for this question are ‘# Define function space’ and ‘# Apply PointSources for lift’