affine_try.py

import pygmsh
import numpy as np
from skfem import *
from skfem.models.elasticity import linear_elasticity, lame_parameters
from skfem.io.meshio import from_meshio
import transformations as tf

BEAM_WIDTH = .635
BEAM_HEIGHT = 0.4318 # .017in=0.4318mm;  .019in=0.4826mm
BEAM_LENGTH = 4

def create_beam():
    geom = pygmsh.opencascade.Geometry(characteristic_length_min=.3, characteristic_length_max=.3 )
    rectangle = geom.add_rectangle([0, 0, 0], BEAM_LENGTH, BEAM_HEIGHT)
    geom.extrude(rectangle, [0, 0, BEAM_WIDTH])

    mesh = pygmsh.generate_mesh(geom)
    return mesh

def affine(x, y, z):
    return np.einsum('ijk->kij', np.einsum('ij,j...', a - np.eye(a.shape[0]), np.stack([x, y, z])) + b)


mesh = create_beam()
m = from_meshio(mesh)

# offset mesh to get origin in the middle of the left side
m.p[1] -= BEAM_HEIGHT/2.0
m.p[2] -= BEAM_WIDTH/2.0

n = m.copy()

# Stiffness Matrix
e1 = ElementTetP1()
e = ElementVectorH1(e1)
ib = InteriorBasis(m, e, MappingIsoparametric(m, e1), 3)
# TODO: Determine Youngs units. Pa? e6 puts it in MPa
K = asm(linear_elasticity(*lame_parameters(4e6, 0.4)), ib)

dofs = {
    'left': ib.get_dofs(lambda x: x[0] == 0.0),
    'right': ib.get_dofs(lambda x: x[0] == BEAM_LENGTH),}

u = np.zeros(K.shape[0])
right_facets = m.facets_satisfying(lambda x: x[0] == BEAM_LENGTH)
right_basis = FacetBasis(m, e, facets=right_facets)
right_dofs = right_basis.get_dofs(right_facets).all()

left_facets = m.facets_satisfying(lambda x: x[0] == 0.0)
left_basis = FacetBasis(m, e, facets=left_facets)
left_dofs = left_basis.get_dofs(left_facets).all()

# Transformation Matrix
# Here I rotate around the z axis (Origin at the center of the left side) but make no translation to approximate
# the functionality of other software packages which only allow translation.
a = tf.rotation_matrix(np.radians(14.32394), [0, 0, 1])[:3,:3]
b = [0,0,0]

# Solve the FEA
u[right_dofs] = L2_projection(affine, right_basis, right_dofs)
I = ib.complement_dofs(dofs)
u = solve(*condense(K, 0 * u, I=I, x=u))

# Displace the Coordinates
sf = 1.0                        # scaling factor
m.p += sf * u[ib.nodal_dofs]    # update the coordinates of the points

# Determine the resultant forces
v = (K @ u)
attempt_1 = np.sum(v[left_dofs].reshape((3, -1)), axis=1)
attempt_2 = np.sum((K @ u).reshape(3, -1), axis=1)


## EX11
v = np.zeros(K.shape[0])
v[dofs['right'].nodal['u^2']] = 1
I = ib.complement_dofs(dofs)
v = solve(*condense(K, 0 * v, I=I, x=v))
sf = 1.0
n.p += sf * v[ib.nodal_dofs]


if __name__ == "__main__":
    print()
    print("***************RESULTS***********************")
    print("Attempt #1:\t\t\t\t\t\t\t {:.2e}, {:.2e}, {:.2e}".format(attempt_1[0], attempt_1[1], attempt_1[2]))
    print("Attempt #2:\t\t\t\t\t\t\t {:.2e}, {:.2e}, {:.2e}".format(attempt_2[0], attempt_2[1], attempt_2[2]))
    print()
    print("SOFTWARE:")
    print("CalculiX:\t\t\t\t\t\t\t {:.2e}, {:.2e}, {:.2e}".format(3.408515e-15, -8.185536e-04, 2.900187e-15))
    print("SolidWorks (Entire Model):\t\t\t {:.2e}, {:.2e}, {:.2e}".format(3.404e-07, -1.9435e-06, 1.7462e-10))
    print("SolidWorks (Left Face):\t\t\t\t {:.2e}, {:.2e}, {:.2e}".format(3.404e-07, 1.2708e-3, 1.7462e-10))
    print("SolidWorks Result (Entire Model):\t {:.2e}".format(1.9731e-06))
    print("SolidWorks Result (Left Face):\t\t {:.2e}".format(1.2708e-3))
    print()
    print("FREE BODY")
    print("SkyCiv: (Beam Calculator)\t\t\t {:.2e} N".format(-.0008))

    from os.path import splitext
    from sys import argv
    from vtkplotter import Plotter
    m.save(splitext(argv[0])[0] + 'try.vtk')
    n.save(splitext(argv[0])[0] + 'ex_11.vtk')


    vp = Plotter(bg=np.array([.223,.223,.223]),axes=2)
    beam1 = vp.load(splitext(argv[0])[0] + 'try.vtk')
    beam1.color('steelblue')
    beam2 = vp.load(splitext(argv[0])[0] + 'ex_11.vtk')
    beam2.color('mediumpurple')
    vp.show([beam1, beam2])