Newer
Older
import numpy as np
import surface_geom_SEM as sgs
from geomdl import exchange
import os
import sys
import cProfile
import pstats
import io
from pstats import SortKey
import subprocess
np.set_printoptions(threshold=sys.maxsize)
Nima
committed
# Global variables
###################################################################
ex = np.array((1, 0, 0), dtype=np.int32)
ey = np.array((0, 1, 0), dtype=np.int32)
ez = np.array((0, 0, 1), dtype=np.int32)
gauss_pw = [[[-0.577350269189626, 1],\
[0.577350269189626, 1]],\
[[-0.774596669241483, 0.555555555555555],\
[0, 0.888888888888888], \
[-0.774596669241483, 0.555555555555555]],\
[[-0.861136311594053, 0.347854845137454],\
[-0.339981043584856, 0.652145154862546],\
[0.339981043584856, 0.652145154862546],\
[0.861136311594053, 0.347854845137454]],
[[-0.906179845938664, 0.236926885056189],\
[-0.538469310105683, 0.478628670499366],\
[0, 0.568888888888889],\
[0.538469310105683, 0.478628670499366],\
[0.906179845938664, 0.236926885056189]]] #to be completed
######################################################################
def profile(func):
'''This function is used for profiling the file.
It will be used as a decorator.
output:
A powershell profling and process.profile file'''
def inner(*args, **kwargs):
pr = cProfile.Profile()
pr.enable()
value = func(*args, **kwargs)
pr.disable()
# s = io.StringIO()
# sortby = SortKey.CUMULATIVE
# ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
# ps.print_stats()
# print(s.getvalue())
pr.dump_stats('process.profile')
return value
return inner
def lbto_pw(data_file_address):
'''This function gets the data from
Mathematica provided .dat file. The data is
lobatto (lbto) points and weights of numeriucal integration
-output:
a 2D np.array in that first column is the
spectral node coordinates and the second column
is the ascociated weights
'''
node_data = np.genfromtxt(data_file_address)
return node_data
def xi_to_uv(xi1, xi2, node_1_u, node_1_v, node_3_u, node_3_v):
'''This function is the mapping form xi1-xi2 which
are parametric space of SEM to u-v space which are
parametric space of the IGA
-output:
u and v according to the input xi1 and xi2 and two oppsite
corners of the rectangle.
'''
if xi1<-1 and xi1>1 and xi2<-1 and xi2>1:
raise ValueError('-1<=xi1<=+1, -1<=xi2<=+1 must be followed')
else:
node_2_u = node_3_u
u = 1/2*((1-xi1)*node_1_u + (1+xi1)*node_2_u)
v = 1/2*((1-xi2)*node_1_v + (1+xi2)*node_3_v)
return u, v
def coorsys_director(surface, u, v):
'''According to P430, finite element
procedures, K. J. Bathe
The output coordinate systems are orhtonormal
ones used to model the rotation of
the director. It is used to define rotations.
-output:
2D (3X3) np.array in that:
the first row is the v1 normalized and
the last row is the director which is already
normalized
'''
vdir = surface.director(u,v)
v1 = (np.cross(ey, vdir))
v1_unit= v1 / np.linalg.norm(v1)
v2_unit = np.cross(vdir, v1_unit)
return np.array([v1_unit, v2_unit, vdir])
def lagfunc(lobatto_pw, xi):
'''Function for generating lagrangian
shape functions
output:
1D np.array in that i-th element
is the value of the i-th lagrangian shape function
at xi. The number of elements in the array obviousluy is
equal to the number of nodes in 1 direction
'''
if xi<-1 and xi>1:
raise ValueError('-1<=xi1<=+1, -1<=xi2<=+1 must be followed')
else:
len = lobatto_pw.shape[0]
# lag_func = np.zeros(len)
# if xi in lobatto_pw[:,0]:
# xi_index = np.where(lobatto_pw[:,0] == xi)
# lag_func[xi_index] = 1
# else: No impresive increase in speed
lag_func = np.ones(len)
for i in range(len):
for j in range(len):
if j != i:
lag_func[i] = lag_func[i] * (xi-lobatto_pw[j, 0]) /\
(lobatto_pw[i, 0]-lobatto_pw[j, 0])
return lag_func
def der_lagfunc_dxi(lobatto_pw, xi):
'''By taking 'ln' functiom from both sides of the definition
of lagrangina shape functions, the derivatives
can be calculated. For example:
https://math.stackexchange.com/questions/809927/first-derivative-of-lagrange-polynomial
however the formulation result in infinity when xi = xj therefore,
it is modified to:
dlag-i-dxi = (sigma (j, j!=i,(PI (k, k!=(i and j), (xi-xk))))/PI(j, j!=i, (xi-xj))
-output:
1D np.array in that i-th element
is the value of the derivatives of i-th lagrangian shape function
w. r. t. xi at xi . The number of element obviousluy is
equal to the number of nodes
'''
len = lobatto_pw.shape[0]
der_lag_dxi = np.zeros(len)
for i in range(len):
numerator_sigma = 0
denominator = 1
for j in range(len):
if j != i:
denominator = denominator * (lobatto_pw[i, 0]-lobatto_pw[j, 0])
numerator_pi = 1
for k in range(len):
if k!=i and k!=j:
numerator_pi = numerator_pi * (xi - lobatto_pw[k, 0])
numerator_sigma = numerator_sigma + numerator_pi
der_lag_dxi [i] = numerator_sigma / denominator
return der_lag_dxi
def der_disp_dxi(surface, lobatto_pw, v1_v2_array, lag_xi1, lag_xi2,\
der_lag_dxi1, der_lag_dxi2, t=0):
'''In this fuctions the derivatives of displacement(ux, uy, uz)
with respect to natural SEM coordinates are calculated. Each node has
5 DOFs: ux_k, uy_k, uz_k, alhpha_k and beta_k
output:
A tuple (der_ux_dxi, der_uy_dxi, der_uz_dxi)in which i-th element
is a (3x3(number of sD-SEM nodes)) np.array.
The below equations show how the tuple elements can be used.
[dux_dxi1, dux_dxi2, dux_dt] = der_ux_dxi .[...ux_k, alpha_k,beta_K...]
[duy_dxi1, duy_dxi2, duy_dt] = der_uy_dxi .[...uy_k, alpha_k,beta_K...]
[duz_dxi1, duz_dxi2, duz_dt] = der_uz_dxi .[...uz_k, alpha_k,beta_K...]
'''
len = lobatto_pw.shape[0]
der_ux_dxi = np.zeros((3, 3*len**2))
der_uy_dxi = np.zeros((3, 3*len**2))
der_uz_dxi = np.zeros((3, 3*len**2))
num = 0
for i in range(len):
for j in range(len):
v1_unit = v1_v2_array[i, j , 0]
v2_unit = v1_v2_array[i, j , 1]
# regarding coorsys_director, v1 is always normal to y axis,
# then the following Eqs. can be simplified, but it is left in
# the general form as a different director coordinate system
# may be used in future.
der_ux_dxi [:,num:num+3] = np.array([[der_lag_dxi1[j] * lag_xi2[i], \
der_lag_dxi1[j] * lag_xi2[i] * t/2 * surface.thk * (-v2_unit[0]),\
der_lag_dxi1[j] * lag_xi2[i] * t/2 * surface.thk * (v1_unit[0])],\
[lag_xi1[j] * der_lag_dxi2[i],\
lag_xi1[j] * der_lag_dxi2[i] * t/2 * surface.thk * (-v2_unit[0]),\
lag_xi1[j] * der_lag_dxi2[i] * t/2 * surface.thk * (v1_unit[0])],\
[0, lag_xi1[j] * lag_xi2[i] * surface.thk/2 * (-v2_unit[0]),\
lag_xi1[j] * lag_xi2[i] * surface.thk/2 * (v1_unit[0])]])
der_uy_dxi [:,num:num+3] = [[der_lag_dxi1[j] * lag_xi2[i], \
der_lag_dxi1[j] * lag_xi2[i] * t/2 * surface.thk * (-v2_unit[1]),\
der_lag_dxi1[j] * lag_xi2[i] * t/2 * surface.thk * (v1_unit[1])],\
[lag_xi1[j] * der_lag_dxi2[i], \
lag_xi1[j] * der_lag_dxi2[i] * t/2 * surface.thk * (-v2_unit[1]),\
lag_xi1[j] * der_lag_dxi2[i] * t/2 * surface.thk * (v1_unit[1])],\
[0, lag_xi1[j] * lag_xi2[i] * surface.thk/2 * (-v2_unit[1]),\
lag_xi1[j] * lag_xi2[i] * surface.thk/2 * (v1_unit[1])]]
Nima
committed
der_uz_dxi [:,num:num+3] = [[der_lag_dxi1[j] * lag_xi2[i], \
der_lag_dxi1[j] * lag_xi2[i] * t/2 * surface.thk * (-v2_unit[2]),\
der_lag_dxi1[j] * lag_xi2[i] * t/2 * surface.thk * (v1_unit[2])],\
[lag_xi1[j] * der_lag_dxi2[i], \
lag_xi1[j] * der_lag_dxi2[i]* t/2*surface.thk*(-v2_unit[2]),\
lag_xi1[j] * der_lag_dxi2[i] * t/2 * surface.thk * (v1_unit[2])],\
[0, lag_xi1[j] * lag_xi2[i] *surface.thk/2 * (-v2_unit[2]),\
lag_xi1[j] * lag_xi2[i] * surface.thk/2 * (v1_unit[2])]]
num = num+3
return (der_ux_dxi, der_uy_dxi, der_uz_dxi)
def der_disp(surface, lobatto_pw, v1_v2_array, lag_xi1, lag_xi2,\
der_lag_dxi1, der_lag_dxi2, node_1_u, node_1_v,\
Nima
committed
node_3_u, node_3_v, jac2, t=0):
'''In this fuctions the derivatives of displacement in x direction
ux, in y direction uy and in z direction uz with respect to
physical coordinates with the help of mapping (jacobian) matrix
are calculated (according to P439, finite element
procedures, K. J. Bathe)
output:
A tuple (der_ux, der_uy, der_uz)in which i-th element
is a 3x3(number of 2D-SEM nodes).
[dux_dx, dux_dy, dux_dz] = der_ux .[...ux_k, alpha_k,beta_K...]
[duy_dx, duy_dy, duy_dz] = der_uy .[...uy_k, alpha_k,beta_K...]
[duz_dx, duz_dy, duz_dz] = der_uz .[...uz_k, alpha_k,beta_K...]'''
len = lobatto_pw.shape[0]
jac_total = np.zeros((3,3))
inv_jac_total = np.zeros((3,3))
der_ux = np.zeros((3, 3*len**2))
der_uy = np.zeros((3, 3*len**2))
der_uz = np.zeros((3, 3*len**2))
Nima
committed
jac1 = np.array([[(node_3_u - node_1_u)/2, 0 , 0],\
[0, (node_3_v - node_1_v)/2, 0],[0, 0, 1]]) #From (xi1, xi2) space to (u,v) space
Nima
committed
# jac2 = surface.ders_uvt(u, v, t)
jac_total = np.matmul(jac1, jac2)
inv_jac_total = np.linalg.inv(jac_total)
der_disp_dxi_var = der_disp_dxi(surface, lobatto_pw, v1_v2_array, lag_xi1, lag_xi2,\
der_lag_dxi1, der_lag_dxi2, t)
der_ux = np.matmul(inv_jac_total, der_disp_dxi_var[0])
der_uy = np.matmul(inv_jac_total, der_disp_dxi_var[1])
der_uz = np.matmul(inv_jac_total, der_disp_dxi_var[2])
return (der_ux, der_uy, der_uz)
def b_matrix(surface, lobatto_pw, v1_v2_array, lag_xi1, lag_xi2,\
der_lag_dxi1, der_lag_dxi2, node_1_u, node_1_v,\
Nima
committed
node_3_u, node_3_v, jac2, t):
'''according to P440, finite element
procedures, K. J. Bathe, B matrix is calculated and
strain = B.d
output:
A 2D numpy array. Its dimension is (6x5(number of nodes))
as each node has 5 DOFs. Strain vector is assumed as
[exx, eyy, ezz, exy, exz, eyz]
'''
len = lobatto_pw.shape[0]
der_disp_var = der_disp(surface, lobatto_pw, v1_v2_array, lag_xi1, lag_xi2,\
der_lag_dxi1, der_lag_dxi2, node_1_u, node_1_v,\
Nima
committed
node_3_u, node_3_v, jac2, t)
der_ux_var = der_disp_var[0]
der_uy_var = der_disp_var[1]
der_uz_var = der_disp_var[2]
b_matrix_var = np.array([[der_ux_var[0,0], 0, 0,
der_ux_var[0,1],der_ux_var[0,2]],
[0, der_uy_var[1,0], 0,
der_uy_var[1,1], der_uy_var[1,2]],
[0, 0, der_uz_var[2,0],
der_uz_var[2,1], der_uz_var[2,2]],
[der_ux_var[1,0], der_uy_var[0,0], 0,
der_ux_var[1,1] + der_uy_var[0,1],
der_ux_var[1,2] + der_uy_var[0,2]],
[0, der_uy_var[2,0], der_uz_var[1,0],
der_uy_var[2,1] + der_uz_var[1,1],
der_uy_var[2,2] + der_uz_var[1,2]],
[der_ux_var[2,0], 0, der_uz_var[0,0],
der_ux_var[2,1] + der_uz_var[0,1],
der_ux_var[2,2]+ der_uz_var[0,2]]])
i = 1
j = 3
len_2 = len**2
while i <= (len_2)-1:
b_local = np.array([[der_ux_var[0,j], 0, 0,
der_ux_var[0,j+1],der_ux_var[0,j+2]],
[0, der_uy_var[1,j], 0,
der_uy_var[1,j+1], der_uy_var[1,j+2]],
[0, 0, der_uz_var[2,j],
der_uz_var[2,j+1], der_uz_var[2,j+2]],
[der_ux_var[1,j], der_uy_var[0,j], 0,
der_ux_var[1,j+1] + der_uy_var[0,j+1],
der_ux_var[1,j+2] + der_uy_var[0,j+2]],
[0, der_uy_var[2,j], der_uz_var[1,j],
der_uy_var[2,j+1] + der_uz_var[1,j+1],
der_uy_var[2,j+2] + der_uz_var[1,j+2]],
[der_ux_var[2,j], 0, der_uz_var[0,j],
der_ux_var[2,j+1] + der_uz_var[0,j+1],
der_ux_var[2,j+2]+ der_uz_var[0,j+2]]])
b_matrix_var = np.append(b_matrix_var, b_local, axis=1)
j += 3
i += 1
return b_matrix_var
Nima
committed
def area(surface, lobatto_pw, node_1_u, node_1_v,\
node_3_u, node_3_v, second_surf_derivatives):
'''This function is only for the test of the
validity of derivatives of mapping and jacobian'''
Nima
committed
det_jac1 =np.linalg.det(np.array([[(node_3_u - node_1_u)/2, 0 , 0],\
[0, (node_3_v - node_1_v)/2, 0],[0, 0, 1]]))
len = lobatto_pw.shape[0]
local_area = 0
i = 0
j = 0
for i in range(len):
Nima
committed
u = xi_to_uv(lobatto_pw[i,0], 0, node_1_u, node_1_v, node_3_u, node_3_v)[0]
for j in range(len):
Nima
committed
v = xi_to_uv( 0, lobatto_pw[j,0], node_1_u, node_1_v, node_3_u, node_3_v)[1]
ders = surface.ders_uvt(0, second_surf_derivatives)
vector_1 = ders[0,:]
vector_2 = ders[1,:]
local_area = local_area +\
np.linalg.norm((np.cross(vector_1,vector_2)))*\
Nima
committed
lobatto_pw[i,1]*lobatto_pw[j,1]
return det_jac1 * local_area
Nima
committed
def physical_crd_xi(surface, xi1, xi2, node_1_u, node_1_v,\
node_3_u, node_3_v, t=0):
u = xi_to_uv(xi1, xi2, node_1_u, node_1_v, node_3_u, node_3_v)[0]
v = xi_to_uv(xi1, xi2, node_1_u, node_1_v, node_3_u, node_3_v)[1]
x = surface.physical_crd(u, v, t)[0]
y = surface.physical_crd(u, v, t)[1]
z = surface.physical_crd(u, v, t)[2]
return x, y, z
def coorsys_material(surface, u, v):
'''According to P438, finite element
procedures, K. J. Bathe Fig. 5.33, an orthonormal coordinate
system, in which the plane stress assumption is valid,
is generated in this function.
-output:
2D (3X3) np.array in that:
the first row is r_hat_unit normalized and
the last row is director which is already
normalized.
surface_der_first = surface.derivatives(u, v, order=1)
g2 = np.array(surface_der_first[0][1])
vdir = surface.director(u,v)
var_1 = np.cross(g2, vdir)
r_hat_unit = var_1 / np.linalg.norm(var_1)
var_2 = np.cross(vdir, r_hat_unit)
s_hat_unit = var_2 / np.linalg.norm(var_2)
return np.array([r_hat_unit, s_hat_unit, vdir])
def mat_transform_matrix(surface, u, v):
'''According to P441, finite element
procedures, K. J. Bathe, a transformation matrix
named in the book as Qsh is generated in this function.
-output:
2D (6X6) np.array
Nima
committed
'''
material_coorsys_vectors = coorsys_material(surface, u, v)
Nima
committed
r_hat_unit = material_coorsys_vectors[0]
s_hat_unit = material_coorsys_vectors[1]
vdir = material_coorsys_vectors[2]
k_1 = np.inner(ex, r_hat_unit) # as both are unit vectors, this is similarity cosine
k_2 = np.inner(ex, s_hat_unit)
k_3 = np.inner(ex, vdir)
m_1 = np.inner(ey, r_hat_unit)
m_2 = np.inner(ey, s_hat_unit)
m_3 = np.inner(ey, vdir)
n_1 = np.inner(ez, r_hat_unit)
n_2 = np.inner(ez, s_hat_unit)
n_3 = np.inner(ez, vdir)
trans_matrix = np.array(((k_1**2, m_1**2, n_1**2, k_1*m_1, m_1*n_1, k_1*n_1),\
(k_2**2, m_2**2, n_2**2, k_2*m_2, m_2*n_2, k_2*n_2),\
(k_3**2, m_3**2, n_3**2, k_3*m_3, m_3*n_3, k_3*n_3),
(2*k_1*k_2, 2*m_1*m_2, 2*n_1*n_2, k_1*m_2 + k_2*m_1,\
m_1*n_2 + m_2*n_1, k_1*n_2 + k_2*n_1),\
(2*k_2*k_3, 2*m_2*m_3, 2*n_2*n_3, k_2*m_3 + k_3*m_2,\
m_2*n_3 + m_3*n_2, k_2*n_3 + k_3*n_2),\
(2*k_1*k_3, 2*m_1*m_3, 2*n_1*n_3, k_1*m_3 + k_3*m_1,\
m_1*n_3 + m_3*n_1, k_1*n_3 + k_3*n_1)))
return trans_matrix
# @profile
def element_stiffness_matrix(surface, lobatto_pw, node_1_u, node_1_v, node_3_u,\
node_3_v, number_gauss_point=2,\
elastic_modulus=3*(10**6), nu=0.3):
'''stiffness matrix and numerical integration
-output:
Element stiffness matrix
'''
len = lobatto_pw.shape[0]
jac1 =np.array(np.array([[(node_3_u - node_1_u)/2, 0 , 0],\
[0, (node_3_v - node_1_v)/2, 0],[0, 0, 1]]))
gauss_pw_nparray = np.array(gauss_pw[number_gauss_point-2])
len_gauss = gauss_pw_nparray.shape[0]
elastic_matrix_planestress = elastic_modulus/(1-nu**2)*\
np.array([[1, nu, 0, 0, 0, 0],[nu, 1, 0, 0, 0, 0],\
[0, 0, 0, 0, 0, 0], [0, 0, 0, (1-nu)/2, 0, 0],
[0, 0, 0, 0, 5/6*(1-nu)/2, 0],[0, 0, 0, 0, 0, 5/6*(1-nu)/2]])
v1_v2_at_lobatto_points = np.zeros((len, len, 2, 3))
for i in range(len):
for j in range(len):
uv_var = xi_to_uv(lobatto_pw[j,0], lobatto_pw[i,0], node_1_u,\
node_1_v, node_3_u, node_3_v)
v1_v2_at_lobatto_points[i, j , 0] = \
coorsys_director(surface, uv_var[0], uv_var[1])[0]
v1_v2_at_lobatto_points[i, j , 1] = \
coorsys_director(surface, uv_var[0], uv_var[1])[1]
stiff_mtx = np.zeros((5*len**2, 5*len**2)) # 5 DOF at each node
for i in range(len):
# print(i, '\n')
xi2 = lobatto_pw[i, 0]
w2 = lobatto_pw[i, 1]
# xi1 = lobatto_pw[i, 0]
# w1 = lobatto_pw[i, 1]
for j in range(len):
xi1 = lobatto_pw[j, 0]
w1 = lobatto_pw[j, 1]
# xi2 = lobatto_pw[j, 0]
# w2 = lobatto_pw[j, 1]
uv_var = xi_to_uv(lobatto_pw[j,0], lobatto_pw[i,0], node_1_u,\
node_1_v, node_3_u, node_3_v)
qsh = mat_transform_matrix(surface, uv_var[0], uv_var[1])
Nima
committed
elastic_matrix = np.transpose(qsh).\
dot(elastic_matrix_planestress).dot(qsh)
lag_xi1 = lagfunc(lobatto_pw, xi1)
lag_xi2 = lagfunc(lobatto_pw, xi2)
der_lag_dxi1 = der_lagfunc_dxi(lobatto_pw, xi1)
der_lag_dxi2 = der_lagfunc_dxi(lobatto_pw, xi2)
second_surf_der = surface.derivatives(uv_var[0], uv_var[1], order=2)
for k in range(len_gauss):
t = gauss_pw_nparray[k, 0]
w3 =gauss_pw_nparray[k, 1]
jac2 = surface.ders_uvt(t, second_surf_der)
b = b_matrix(surface, lobatto_pw, v1_v2_at_lobatto_points,\
lag_xi1, lag_xi2,der_lag_dxi1,\
der_lag_dxi2, node_1_u, node_1_v,\
Nima
committed
node_3_u, node_3_v, jac2, t)
# btr_d_b = np.transpose(b).dot(elastic_matrix).dot(b)
btr_d_b = np.transpose(b) @ elastic_matrix @ b
stiff_mtx = stiff_mtx + np.linalg.det(jac1)*np.linalg.det(jac2)\
*btr_d_b*w1*w2*w3 #1/4 for (xi1, xi2) space to (u,v) space
def check_symmetric(a, rtol=1e-05, atol=1e-05):
return np.allclose(a, a.T, rtol=rtol, atol=atol)
def boundary_condition(lobatto_pw):
# bc_h_bott = [0, 0, 1, 1, 1] #cantilever Plate.zero means clamped DOF
# bc_h_top = [0, 0, 0, 0, 0]
# bc_v_left = [0, 0, 1, 1, 1]
# bc_v_right = [0, 0, 1, 1, 1]
# bc_h_bott = [1, 1, 1, 1, 1] #cantilever quarter cylinder.zero means clamped DOF
# bc_h_top = [0, 0, 0, 0, 0]
# bc_v_left = [1, 1, 1, 1, 1]
# bc_v_right = [1, 1, 1, 1, 1]
bc_h_bott = [1, 0, 1, 0, 1] #pinched shell. zero means clamped DOF
bc_h_top = [0, 1, 0, 1, 0]
bc_v_left = [1, 1, 0, 1, 0] #xi1 = 1, director direction is inwards
bc_v_right = [0, 1, 1, 1, 0]
# bc_h_bott = [0, 0, 0, 0, 0] #Total clamped. zero means clamped DOF
# bc_h_top = [0, 0, 0, 0, 0]
# bc_v_left = [0, 0, 0, 0, 0]
# bc_v_right = [0, 0, 0, 0, 0]
# bc_h_bott = [0, 0, 0, 0, 0] # Finding stiffness zero means clamped DOF
# bc_h_top = [0, 0, 0, 0, 0]
# bc_v_left = [0, 0, 1, 0, 0]
# bc_v_right = [0, 0, 0, 0, 0]
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
len = lobatto_pw.shape[0]
bc_line_h_bott = np.zeros(5*len) #bottom horizontal line
bc_line_h_top = np.zeros(5*len) #top horizontal line
bc_line_v_left = np.zeros(5*(len - 2))
bc_line_v_right = np.zeros(5*(len - 2))
i = 1
c_node = len * (len-1)
while i <= len:
bc_line_h_bott[5*i - 5] = (1 - bc_h_bott[0]) * (5*i-4)
bc_line_h_bott[5*i - 4] = (1 - bc_h_bott[1]) * (5*i-3)
bc_line_h_bott[5*i - 3] = (1 - bc_h_bott[2]) * (5*i-2)
bc_line_h_bott[5*i - 2] = (1 - bc_h_bott[3]) * (5*i-1)
bc_line_h_bott[5*i - 1] = (1 - bc_h_bott[4]) * (5*i)
bc_line_h_top[5*i - 5] = (1 - bc_h_top[0]) * (5*(c_node + i) - 4)
bc_line_h_top[5*i - 4] = (1 - bc_h_top[1]) * (5*(c_node + i) - 3)
bc_line_h_top[5*i - 3] = (1 - bc_h_top[2]) * (5*(c_node + i) - 2)
bc_line_h_top[5*i - 2] = (1 - bc_h_top[3]) * (5*(c_node + i) - 1)
bc_line_h_top[5*i - 1] = (1 - bc_h_top[4]) * (5*(c_node + i))
i += 1
i = 1
while i <= len-2:
bc_line_v_left[5*i - 5] = (1 - bc_v_left[0]) * (5*(i*len + 1) - 4)
bc_line_v_left[5*i - 4] = (1 - bc_v_left[1]) * (5*(i*len + 1) - 3)
bc_line_v_left[5*i - 3] = (1 - bc_v_left[2]) * (5*(i*len + 1) - 2)
bc_line_v_left[5*i - 2] = (1 - bc_v_left[3]) * (5*(i*len + 1) - 1)
bc_line_v_left[5*i - 1] = (1 - bc_v_left[4]) * (5*(i*len + 1))
bc_line_v_right[5*i - 5] = (1 - bc_v_right[0]) * (5*(i*len + len) - 4)
bc_line_v_right[5*i - 4] = (1 - bc_v_right[1]) * (5*(i*len + len) - 3)
bc_line_v_right[5*i - 3] = (1 - bc_v_right[2]) * (5*(i*len + len) - 2)
bc_line_v_right[5*i - 2] = (1 - bc_v_right[3]) * (5*(i*len + len) - 1)
bc_line_v_right[5*i - 1] = (1 - bc_v_right[4]) * (5*(i*len + len))
i += 1
bc_1 = np.concatenate((bc_line_h_bott, bc_line_h_top, bc_line_v_left,\
bc_line_v_right))
#####################################################
'''Only for pinched shell'''
node_3_number = (len-1)*len + 1 #xi1 = -1, xi2 = 1
node_4_number = node_3_number + len - 1 #xi1 = 1, xi2 = 1
bc_node_1 = np.array([0, 2, 3, 4, 5])
bc_node_2 = np.array([5*len-4, 5*len-3, 0, 5*len-1, 5*len])
bc_node_3 = np.array([5*node_3_number-4, 0, 5*node_3_number-2, 0, 5*node_3_number])
bc_node_4 = np.array([5*node_4_number-4, 0, 5*node_4_number-2, 0, 5*node_4_number])
for i in range(5):
bc_1[i] = bc_node_1[i]
bc_1[5*len - 5 + i] = bc_node_2[i]
bc_1[5*(len + 1) - 5 + i] = bc_node_3[i]
bc_1[5*(len + 1 + len - 1) - 5 + i] = bc_node_4[i]
#####################################################
bc_2 = np.sort(np.delete(bc_1, np.where(bc_1 == 0))) - 1 # in python arrays start from zero
return (bc_2.astype(int))
Nima
committed
Nima
committed
def stiffness_matrix_bc_applied(stiff_mtx, bc):
stiff = np.delete(np.delete(stiff_mtx, bc, 0), bc, 1)
return stiff
if __name__ == '__main__':
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
############## Visualization
# surf_cont = multi.SurfaceContainer(data)
# surf_cont.sample_size = 30
# surf_cont.vis = vis.VisSurface(ctrlpts=False, trims=False)
# surf_cont.render()
########################## TESTs
# Test for area and mapping matrix
# area1 = area(surfs, lobatto_pw)
# print("area is :", area1)
#Test for lag_func
# lobatto_pw = lbto_pw("node_weight.dat")
# lag = lagfunc(lobatto_pw, -0.8)
# tt =2
# Test for coorsys_material
# coor = coorsys_material(surfs, 0.5, -0.3) #only the dimensions approved
# print(coor, '\n')
#Test for coorsys_material
# coor = mat_transform_matrix(surfs, 0.5, -0.2)
# print(coor, '\n')
# print("area is :", area1)
# print(physical_crd_xi(surfs, 1, 1))
# print(physical_crd_xi(surfs, -1, 1))
# print(physical_crd_xi(surfs, -1, -1))
# print(physical_crd_xi(surfs, 1, -1))
# Test for tangent and director vector
# len = lobatto_pw.shape[0]
# for i in range(len):
# for j in range(len):
# xi1 = lobatto_pw[j, 0]
# xi2 = lobatto_pw[i, 0]
# xyz = physical_crd_xi(surfs, xi1, xi2)
# vectors = coorsys_material(surfs, xi1, xi2)
# print('x, y, z= ', xyz)
# print('vectors = ', vectors, '\n')
#Test for xi_to_uv
# print(xi_to_uv(-1, 0.5, 0, 0, 5, 5))
#########################################################
data = exchange.import_json("pinched_shell.json") # pinched_shell.json rectangle_cantilever square square_kninsertion generic_shell_kninsertion foursided_curved_kninsertion foursided_curved_kninsertion2 rectangle_kninsertion
surfs = sgs.SurfaceGeo(data, 0, 3)
lobatto_pw = lbto_pw("node_weight.dat")
node_1_u = 0
node_1_v = 0
node_3_u = 1
node_3_v = 1
time1 = time.time()
stiff = element_stiffness_matrix(surfs, lobatto_pw, node_1_u, node_1_v, node_3_u,\
node_3_v)
bc = boundary_condition(lobatto_pw)
print(bc)
stiff_bc = stiffness_matrix_bc_applied(stiff, bc)
# len = lobatto_pw.shape[0]
# print('number of nodes:',len,'\n')
load = np.zeros(np.shape(stiff)[0])
p =1000
# # # load[5*int((len-1)/2) + 3 - 1] = p #for cantilever
# # load[8:(5*int((len))-1 - 1):5] = p #for cantilever subjected to moment at free end
# # load[3] = p/2 #for cantilever subjected to moment at free end
# # load[(5*int((len))-1 - 1)] =p/2 #for cantilever subjected to moment at free end
# # # load[5*int((len-1)/2*len+(len-1)/2) + 3 - 1] = 1
load[0] = p #for pinched
load = np.delete(load, bc, 0)
d = np.linalg.solve(stiff_bc, load)
# # print(d)
# # print("displacement:", d[int(3*(len-1)/2)+1-1], '\n') #for cantilever subjected to moment at free end
# time2 = time.time()
# print( "time is:",time2-time1, '\n')
print("displacement:", d[np.where(load==p)[0][0]], '\n')
# print(f'condition number of the stiffness matrix : {np.linalg.cond(stiff_bc)}\n')
# # print("determinant of stiff matrix is", np.linalg.det(stiff_bc))
# print("End")
subprocess.call("C:\\Users\\azizi\Desktop\\DFG_Python\\.P394\\Scripts\\snakeviz.exe process.profile ", \
shell=False)
# "editor.hover.enabled": false,
# "editor.hover.sticky": false,
# "editor.hover.above": false,
# "editor.parameterHints.enabled": false,
# "editor.semanticHighlighting.enabled": true,