Skip to content
Snippets Groups Projects
curvature_calculation.py 2.3 KiB
Newer Older
from geomdl import exchange
import surface_geom_SEM as sgs
import numpy as np
import element_stiff_matrix as esm
###########################################################################################

with open("crv_intpoint.dat", "w") as file1:
    pass
data = exchange.import_json("double-curved-free-form_5B.json") #double-curved-free-form.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.1
surfs = sgs.SurfaceGeo(data, 0, thk)
lobatto_pw_all =esm.lbto_pw("node_weight_all.dat")
i_main = 14
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, :]
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 = 0 # t =xi3
# area1= esm.area(surfs, lobatto_pw, node_1_u, node_1_v,\
#                 node_3_u, node_3_v, t)
# print("area is :", area1)
crv_matrix_intpoint = np.zeros(5)
# crv_matrix_intpoint = np.zeros(3)
u_vector =np.zeros(dim)
# for i in range(dim):
print(surfs.derivatives(0 ,1 , order=2)[1][0]) 
print(surfs.derivatives(0 ,1 , order=2)[0][1])
xi2 = -1 # lobatto_pw[i, 0]
v = esm.xi_to_uv(0, xi2, node_1_u, node_1_v, node_3_u, node_3_v)[1]
lag_xi2 = esm.lagfunc(lobatto_pw, xi2)
der_lag_dxi2 = esm.der_lagfunc_dxi(lobatto_pw, xi2)
for j in range(dim):
    xi1 = lobatto_pw[j, 0]
    u = esm.xi_to_uv(xi1, 0, node_1_u, node_1_v, node_3_u, node_3_v)[0]
    lag_xi1 = esm.lagfunc(lobatto_pw, xi1)
    der_lag_dxi1 = esm.der_lagfunc_dxi(lobatto_pw, xi1)
    second_surf_der = surfs.derivatives(u ,v , order=2)
    crv_gauss = surfs.curvature_mtx(second_surf_der)
    crv_matrix_intpoint[0] = xi1
    crv_matrix_intpoint[1] = xi2
    crv_matrix_intpoint[2] = u
    crv_matrix_intpoint[3] = v
    crv_matrix_intpoint[4] = -crv_gauss
    u_vector[j] = u
    with open("crv_intpoint.dat", "a") as file1:
        np.savetxt(file1, np.array([crv_matrix_intpoint]))
print(u_vector)