Newer
Older
import numpy as np
import surface_geom_SEM as sgs
import surface_nurbs_isoprm as snip
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)
if abs(abs(vdir[1])-1) > 0.00001 :
v1 = (np.cross(ey, vdir))
v1_unit= v1 / np.linalg.norm(v1)
v2_unit = np.cross(vdir, v1_unit)
else:
if vdir[1] > 0:
v1_unit = np.array((1, 0, 0))
v2_unit = np.array((0, 0, -1))
vdir = np.array((0, 1, 0))
else:
v1_unit = np.array((1, 0, 0))
v2_unit = np.array((0, 0, 1))
vdir = np.array((0, 1, 0))
# vdir = surface.director(u,v) # Lamina coordinate system
# surface_der = surface.derivatives(u, v, 1)
# g1 = np.array(surface_der[1][0])
# g2 = np.array(surface_der[0][1])
# g1_unit = g1 / np.linalg.norm(g1)
# g2_unit = g2 / np.linalg.norm(g2)
# A_xi1 = (0.5*(g1_unit + g2_unit))/np.linalg.norm(0.5*(g1_unit + g2_unit))
# v_var_1 = np.cross(vdir, A_xi1)
# A_xi2 = (v_var_1) / np.linalg.norm(v_var_1)
# v1_unit = (2**0.5)/2 * (A_xi1 - A_xi2)
# v2_unit = (2**0.5)/2 * (A_xi1 + A_xi2)
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 SEM nodes^2)) 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] = [[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,\
node_3_u, node_3_v, jac2, t):
'''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))
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 np.identity(3)
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, eyz, exz]
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, t):
'''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]
second_surf_der = surfs.derivatives(u, v, 2)
ders = surface.ders_uvt(t, second_surf_der)
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
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
def volum(surface, lobatto_pw, number_gauss_point, node_1_u, node_1_v,\
node_3_u, node_3_v):
'''This function is only for the test of the
validity of derivatives of mapping and jacobian'''
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]
gauss_pw_nparray = np.array(gauss_pw[number_gauss_point-2])
len_gauss = gauss_pw_nparray.shape[0]
# i = 0
# j = 0
# k = 0
vol = 0
for i in range(len):
xi2 = lobatto_pw[i, 0]
w2 = lobatto_pw[i, 1]
v = xi_to_uv(0, lobatto_pw[i,0], node_1_u, node_1_v, node_3_u, node_3_v)[1]
# print('v is: ',v)
for j in range(len):
xi1 = lobatto_pw[j, 0]
w1 = lobatto_pw[j, 1]
u = xi_to_uv(lobatto_pw[j,0], 0, node_1_u, node_1_v, node_3_u, node_3_v)[0]
# print('u is: ', u)
second_surf_der = surfs.derivatives(u, v, 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)
vol = vol + np.linalg.det(jac2)* det_jac1 * w1 * w2 *w3
return vol
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]
physical_coor_var = surface.physical_crd(u, v, t)
x = physical_coor_var[0]
y = physical_coor_var[1]
z = physical_coor_var[2]
# 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, elastic_modulus=3*(10**6), \
nu=0.3, number_gauss_point=2):
'''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]])) #np.identity(3)
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
# with open ('jac_semn.dat', 'a') as jac_file:
# np.savetxt(jac_file, (node_1_u, node_1_v, node_3_u,\
# node_3_v))
# jac_file.write("\n")
# with open("b_semn.dat", 'w') as b_file:
# pass
# with open ('qsh_semn.dat', 'w') as qsh_file:
# pass
# with open("kwojac_semn.dat", 'w') as k_file:
# pass
# with open("director_unit_semn.dat", 'w') as du:
# pass
# with open("k_semn.dat", 'w') as k_file:
# pass
# with open("jac_semn.dat", 'w') as jac_file:
# pass
surface_isoprm = snip.SurfaceGeneratedSEM(surface, lobatto_pw, node_1_u,\
node_1_v, node_3_u, node_3_v)
# coorsys_tanvec_mtx = surface_isoprm.coorsys_director_tanvec_allnodes()
for i in range(len):
# print(i, '\n')
xi2 = lobatto_pw[i, 0]
w2 = lobatto_pw[i, 1]
lag_xi2 = lagfunc(lobatto_pw, xi2)
der_lag_dxi2 = der_lagfunc_dxi(lobatto_pw, xi2)
# 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])
# with open ('qsh_semn.dat', 'a') as qsh_file:
# np.savetxt(qsh_file, qsh)
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)
crv_content = surface.curvature_mtx(second_surf_der)
gauss_crv = crv_content
with open ('gauss_crv_semn.dat', 'a') as g_file:
np.savetxt(g_file, [gauss_crv])
# g_file.write("\n")
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) #np.identity(3)
# jac2 = surface.jacobian_semi(coorsys_tanvec_mtx, lobatto_pw, i, j, \
# t, second_surf_der, lag_xi1, lag_xi2, der_lag_dxi1, der_lag_dxi2)
################
with open ('jac_semn.dat', 'a') as jac_file:
np.savetxt(jac_file, [np.linalg.det(jac1)*np.linalg.det(jac2)])
jac_file.write("\n")
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)
# if i == 2 and j == 1:#UV
# print(f"[{i}, {j}]= ", b[0, 45])
# print(f"[{i}, {j}]= ", b[0, 46])
# print(f"[{i}, {j}]= ", b[0, 47])
# print(f"[{i}, {j}]= ", b[0, 48])
# print(f"[{i}, {j}]= ", b[0, 49])
# print(f"ey[{i}, {j}]= ", b[1, 45])
# if i == 1 and j == 2:#UV changed
# print(f"[{i}, {j}]= ", b[0, 30])
# print(f"[{i}, {j}]= ", b[0, 31])
# print(f"[{i}, {j}]= ", b[0, 32])
# print(f"[{i}, {j}]= ", b[0, 33])
# print(f"[{i}, {j}]= ", b[0, 34])
# print(f"ey[{i}, {j}]= ", b[1, 31])
# btr_d_b = np.transpose(b).dot(elastic_matrix).dot(b)
btr_d_b = np.transpose(b) @ elastic_matrix @ b
# with open ('b_semn.dat', 'a') as b_file:
# np.savetxt(b_file, btr_d_b)
stiff_mtx = stiff_mtx + np.linalg.det(jac1)*np.linalg.det(jac2)\
Nima
committed
*btr_d_b*w1*w2*w3
# with open ('kwojac_semn.dat', 'a') as k_file:
# np.savetxt(k_file, stiff_mtx)
# with open ('k_semn.dat', 'a') as k_file:
# np.savetxt(k_file, stiff_mtx)
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]
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
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__':
############## 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
data = exchange.import_json("double-curved-free-form_5B.json") # curved_beam_lineload_2_kninsertion curved_beam_lineload_2 pinched_shell_kninsertion_changedeg.json pinched_shell.json rectangle_cantilever square square_kninsertion generic_shell_kninsertion foursided_curved_kninsertion foursided_curved_kninsertion2 rectangle_kninsertion
# visualization(data)
thk = 0.5
surfs = sgs.SurfaceGeo(data, 0, thk)
lobatto_pw_all =lbto_pw("node_weight_all.dat")
i_main = 8
if i_main == 1:
lobatto_pw = lobatto_pw_all[1:3,:]
else:
index = np.argwhere(lobatto_pw_all==i_main)
lobatto_pw = lobatto_pw_all[index[0, 0] + 1:\
index[0, 0] + (i_main+1) +1, :]
Nima
committed
dim = lobatto_pw.shape[0]
# jac_lobatto_node = np.zeros((dim, dim, 3, 3))
node_1_u = 0
node_1_v = 0
node_3_u = 1
node_3_v = 1
t = 1
Nima
committed
with open("jacobian.dat", "w") as f:
pass
second_surf_der = surfs.derivatives(0.5, 0.5, order=2)
# print(second_surf_der)\
print('g1 = ', second_surf_der[1][0])
print('\ng2 = ', second_surf_der[0][1],'\n\n')
print('d2phi/du^2 = ', second_surf_der[2][0])
print('\nd2phi/duv = ', second_surf_der[1][1])
print('\nd2phi/dv^2 = ', second_surf_der[0][2])
Nima
committed
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
for i in range(dim):
for j in range(dim):
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)
second_surf_der = surfs.derivatives(uv_var[0], uv_var[1], order=2)
g1 = np.array(second_surf_der[1][0])
g2 = np.array(second_surf_der[0][1])
dg1_du = np.array((second_surf_der[2][0][0],\
second_surf_der[2][0][1], second_surf_der[2][0][2]))
dg1_dv = np.array((second_surf_der[1][1][0],\
second_surf_der[1][1][1], second_surf_der[1][1][2]))
dg2_dv = np.array((second_surf_der[0][2][0], \
second_surf_der[0][2][1], second_surf_der[0][2][2]))
dg2_du = dg1_dv
vdir_notunit = np.cross(g1, g2)
dvdir_notunit_du = np.cross(dg1_du, g2) + np.cross(g1, dg2_du)
norm_vdir_notunit = np.linalg.norm(vdir_notunit)
part_1 = np.array([[dg1_du[0], g2[0], vdir_notunit[0]],
[dg1_du[1], g2[1], vdir_notunit[1]],
[dg1_du[2], g2[2], vdir_notunit[2]]])
part_2 = np.array([[g1[0], dg2_du[0], vdir_notunit[0]],
[g1[1], dg2_du[1], vdir_notunit[1]],
[g1[2], dg2_du[2], vdir_notunit[2]]])
part_3 = np.array([[g1[0], g2[0], dvdir_notunit_du[0]],
[g1[1], g2[1], dvdir_notunit_du[1]],
[g1[2], g2[2], dvdir_notunit_du[2]]])
dnorm_dvdir_notunit_du = (np.linalg.det(part_1) +\
np.linalg.det(part_2) +\
np.linalg.det(part_3))/\
(2*norm_vdir_notunit)
dvdir_unit_du = (dvdir_notunit_du*norm_vdir_notunit -\
vdir_notunit*dnorm_dvdir_notunit_du)/norm_vdir_notunit**2
dvdir_notunit_dv = np.cross(dg1_dv, g2) + np.cross(g1, dg2_dv)
part_1 = np.array([[dg1_dv[0], g2[0], vdir_notunit[0]],
[dg1_dv[1], g2[1], vdir_notunit[1]],
[dg1_dv[2], g2[2], vdir_notunit[2]]])
part_2 = np.array([[g1[0], dg2_dv[0], vdir_notunit[0]],
[g1[1], dg2_dv[1], vdir_notunit[1]],
[g1[2], dg2_dv[2], vdir_notunit[2]]])
part_3 = np.array([[g1[0], g2[0], dvdir_notunit_dv[0]],
[g1[1], g2[1], dvdir_notunit_dv[1]],
[g1[2], g2[2], dvdir_notunit_dv[2]]])
dnorm_dvdir_notunit_dv = (np.linalg.det(part_1) +\
np.linalg.det(part_2) +\
np.linalg.det(part_3))/\
(2*norm_vdir_notunit)
dvdir_unit_dv = (dvdir_notunit_dv*norm_vdir_notunit -\
vdir_notunit*dnorm_dvdir_notunit_dv)/ norm_vdir_notunit**2
vdir_unit = vdir_notunit / np.linalg.norm(vdir_notunit)
# dx_du = g1[0] + t/2 * thk * dvdir_unit_du[0]
# dy_du = g1[1] + t/2 * thk * dvdir_unit_du[1]
# dz_du = g1[2] + t/2 * thk * dvdir_unit_du[2]
# dx_dv = g2[0] + t/2 * thk * dvdir_unit_dv[0]
# dy_dv = g2[1] + t/2 * thk * dvdir_unit_dv[1]
# dz_dv = g2[2] + t/2 * thk * dvdir_unit_dv[2]
# dx_dt = 1/2 * thk * vdir_unit[0]
# dy_dt = 1/2 * thk * vdir_unit[1]
# dz_dt = 1/2 * thk * vdir_unit[2]
dx_du = dvdir_unit_du[0]
dy_du = dvdir_unit_du[1]
dz_du = dvdir_unit_du[2]
Nima
committed
dx_dv = dvdir_unit_dv[0]
dy_dv = dvdir_unit_dv[1]
dz_dv = dvdir_unit_dv[2]
Nima
committed
dx_dt = vdir_unit[0]
dy_dt = vdir_unit[1]
dz_dt = vdir_unit[2]
Nima
committed
ders = np.array(((dx_du, dy_du, dz_du),(dx_dv, dy_dv, dz_dv),\
(dx_dt, dy_dt, dz_dt)))
# jac_lobatto_node[i, j, : , : ] = ders
jac1 =np.array(np.array([[1/2, 0 , 0],\
[0, 1/2, 0],[0, 0, 1]]))
ders_total = jac1 @ ders
with open ('jacobian.dat', 'a') as jac_file:
np.savetxt(jac_file, ders_total)
jac_file.write("\n")
Nima
committed
#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 lag_func
lobatto_pw = lbto_pw("node_weight.dat")
der_lag = der_lagfunc_dxi(lobatto_pw, -1 )
# print(der_lag, '\n')
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
#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))
#########################################################
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
# data = exchange.import_json("sphere_clean_notrim.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)
# a0 = coorsys_director(surfs, 0.0 , 0.0)
# a1 = coorsys_director(surfs, 1.0 , 0.0)
# a2 = coorsys_director(surfs, 1.0 , 1.0)
# a3 = coorsys_director(surfs, 0.0 , 1.0)
# print(a0, '\n', a1, '\n', a2, '\n',a3)
# lobatto_pw = lbto_pw("node_weight.dat")
# print(area(surfs, lobatto_pw, 0, 0,\
# 1, 1))
# 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,