Beam Reactions

#!/usr/bin/env python
""" Cantilever Timoshenko beam problem solved using
continuous-displacement/constant-shear finite element
method. Cantilever type boundary conditions with point
load applied.

Author: Jack S. Hale 2014
Email: [email protected]

FE discretisation: Jack S. Hale PhD p. 114
Beam problem: Jack S. Hale PhD p. 124
Analytical solution: Jack S. Hale PhD p.125
"""

import numpy as np
from dolfin import *
# from dolfin_error_norms import relative_error, l2, h1


parameters['form_compiler']['representation'] = 'quadrature'
parameters['form_compiler']['quadrature_degree'] = 2



def main(nx=24, degree=2, epsilon=0.0001):
    mesh = UnitIntervalMesh(nx)

    R = FunctionSpace(mesh, "CG", degree)
    V_3 = FunctionSpace(mesh, "CG", degree)
    S = FunctionSpace(mesh, "DG", degree - 1)

    # the method works with no locking
    L = 1.0
    p_tilde = 3.0

    # Inverse of shear constant
    C_m = Constant(epsilon**2.0)
    # Bending constants
    C_b = Constant(L**2.0)

    # Rotation, Displacement, Shear Stress
    U = MixedFunctionSpace([R, V_3, S])
    R, V_3, S = U.split()

    # Rotation, Displacement, Shear Stress test functions
    theta, z_3, gamma = TrialFunctions(U)
    # Rotation, Displacement, Shear Stress trial functions
    eta, y_3, psi = TestFunctions(U)

    # bilinear form Jack S. Hale PhD p.110
    A = C_b*inner(grad(theta), grad(eta))*dx  + \
        inner(gamma, y_3.dx(0) - eta)*dx + \
        inner(z_3.dx(0) - theta, psi)*dx - \
        C_m*inner(gamma, psi)*dx

    # linear form
    L = Constant(0.0)*y_3*dx

    # Fixed (built-in) support on left
    def left_boundary(x, on_boundary):
        tol = 1E-14
        return on_boundary and np.abs(x[0]) < tol

    zero = Constant(0.0)
    bc1 = DirichletBC(V_3, zero, left_boundary)
    bc2 = DirichletBC(R, zero, left_boundary)
    bcs = [bc1, bc2]

    # point load on right
    end_point = Point(1.0)
    f = PointSource(V_3, end_point, p_tilde)

    # Not sure if it is possible to use the PointSource class
    # and LinearVariationalProblem so assemble system manually
    u_h = Function(U)
    A_matrix = assemble(A)
    b_vector = assemble(L)
    for bc in bcs:
        bc.apply(A_matrix)
        bc.apply(b_vector)
    # Apply point source to rhs vector
    f.apply(b_vector)

    # solve
    solver = LUSolver(A_matrix)
    solver.solve(u_h.vector(), b_vector)

    theta_h, z_3h, gamma_h = u_h.split()

    # Jack S. Hale PhD p.125
    z_3 = Expression('0.5*(2.0 - 3.0*(1.0 - x[0]) + pow((1.0 - x[0]), 3) + 6.0*pow(epsilon,2)*x[0])',
                     epsilon=epsilon)
    theta = Expression('1.5*(1.0 - pow((1.0 - x[0]),2))')

    V_e = FunctionSpace(mesh, "CG", 3)

    results = {}
    results['z3_tip_error'] = z_3h(1.0)/z_3(1.0)
    results['z3_tip'] = z_3h(1.0)
    results['theta_tip_error'] = theta_h(1.0)/theta(1.0)
    # results['z3_l2'] = relative_error(z_3, z_3h, V_e, norm=l2)
    # results['z3_h1'] = relative_error(z_3, z_3h, V_e, norm=h1)
    # results['theta_l2'] = relative_error(theta, theta_h, V_e, norm=l2)
    # results['theta_h1'] = relative_error(theta, theta_h, V_e, norm=h1)
    results['z_3h'] = z_3h
    results['theta_h'] = theta
    results['gamma_h'] = gamma_h
    results['dofs'] = U.dim()

    reaction = Function(U)
    reaction.vector()[:] = A_matrix * u_h.vector()
    reaction_forces = reaction.split(deepcopy=True)[1]
    plot(reaction_forces); interactive()
    matrix = assemble(A)
    reaction.vector()[:] = matrix * u_h.vector()
    reaction_forces = reaction.split(deepcopy=True)[1]
    plot(reaction_forces); interactive()
    return results


if __name__ == "__main__":
    results = main(nx=20, degree=2)
    print results

    # Uncomment these lines to plot results
    # plot(results['z_3h'])
    # interactive()