import numpy as np
import os
from geomdl import exchange
import time
import surface_geom_SEM as sgs
import element_stiff_matrix as esm
import global_stiff_matrix as gsm
import global_load_vector_uniform as glv
import gnuplot_automated as afgnu
import global_load_uni_lineload as glineloadv
#INPUT *************************************************************
u_analytic = 0.942
elastic_modulus = 1000
nu = 0.
thk =0.001
bc_h_bott = [0, 0, 0, 0, 0] # zero means clamped DOF
bc_h_top = [1, 1, 1, 1, 1]
bc_v_left = [1, 1, 1, 1, 1]
bc_v_right = [1, 1, 1, 1, 1]
loaded_side = ['top'] # line load is applied in which side in uv space
left_lineload = [0.1*thk**3, 0, 0] #load values in global x, y, z directions
right_lineload = [0.1*thk**3, 0, 0]
bottom_lineload = [0.1*thk**3, 0, 0]
top_lineload = [0.1*thk**3, 0, 0]
#*************************************************************
# os.system('cls')
print("\nImporting Lobatto points and weights from data-base ...")
# time.sleep(1)
lobatto_pw_all = esm.lbto_pw("node_weight_all.dat")
print("\nImporting geometry from json file ...")
# time.sleep(1)
data = exchange.import_json("curved_beam_lineload_2.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)
surfs = sgs.SurfaceGeo(data, 0, thk)
# p_1 = surfs.physical_crd(0., 0.)
# p_2 = surfs.physical_crd(1., 0.)
# p_3 = surfs.physical_crd(1., 1.)
# print("p_1:", p_1, " p_2:", p_2, " p_3:", p_3)
print("\nSelecting the type of refinement: \nPlease enter 'p' for p_refinement and 'h' for h_refinement")
refine_type = input()
# p_1 = surfs.physical_crd(0., 0.)
# p_2 = surfs.physical_crd(1., 0.)
# p_3 = surfs.physical_crd(1., 1.)
# print("p_1:", p_1, " p_2:", p_2, " p_3:", p_3)
if refine_type =="p":
min_order_elem = int(input("\nEnter the minimum order of elements (max order = 30):\n"))
max_order_elem = int(input("Enter the maximum order of elements:\n"))
num_elem_u_auto = int(input("\nEnter the number of elements in u direction:\n"))
num_elem_v_auto = int(input("Enter the number of elements in v direction:\n"))
print("\nEnter the order of continuity at knots to be used for auto detection of elements boundaries in u direction ...")
print("The default value is '1'")
c_order_u =int(input())
print("Enter the order of continuity at knots to be used for auto detection of elements boundaries in v direction ...")
print("The default value is '1'")
c_order_v =int(input())
u_manual = np.linspace(0, 1, num_elem_u_auto + 1) #np.linspace(a, b, c) divide line bc to c-1 parts or add c points to it.
v_manual = np.linspace(0, 1, num_elem_v_auto + 1)
print("\nProgram starts to autodetect the element boundaries and generate mesh ...")
mesh = gsm.mesh_func(surfs, u_manual, v_manual, c_order_u, c_order_v)
element_boundaries_u = mesh[0]
element_boundaries_v = mesh[1]
order_displm_array = np.zeros((max_order_elem - min_order_elem + 1, 2))
time_assembling = np.zeros((max_order_elem - min_order_elem + 1, 2))
time_solver = np.zeros((max_order_elem - min_order_elem + 1, 2))
dof_displm_array = np.zeros((max_order_elem - min_order_elem + 1, 2)) ######################################
dof_time_assembling = np.zeros((max_order_elem - min_order_elem + 1, 2)) ######################################
dof_time_solver = np.zeros((max_order_elem - min_order_elem + 1, 2)) ######################################
cond_elem =np.zeros((max_order_elem - min_order_elem + 1, 2))
order_counter = 0
i_main = min_order_elem
while i_main <= max_order_elem:
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, :]
bc = gsm.global_boundary_condition(lobatto_pw, bc_h_bott, bc_h_top,\
bc_v_left, bc_v_right, element_boundaries_u,\
element_boundaries_v)
print("\n\n\nOrder of elements: {}".format(str(i_main)))
print("\nAssembling global stiffness matrix ...")
t_1_assembling = time.perf_counter()
k_global = gsm.global_stiffness_matrix(surfs, lobatto_pw, \
element_boundaries_u, element_boundaries_v, elastic_modulus, nu)
t_2_assembling = time.perf_counter()
k_global_bc = esm.stiffness_matrix_bc_applied(k_global, bc)
print("\nAssembling global load vector ...")
global_load = glineloadv.global_lineload_vector(surfs, lobatto_pw, \
loaded_side, left_lineload, right_lineload, \
bottom_lineload, top_lineload, element_boundaries_u, \
element_boundaries_v)
global_load_bc = np.delete(global_load, bc, 0)
# print('bc:', bc)
# print(global_load_bc)
print("\nLinear solver in action! ...")
t_1_solver = time.perf_counter()
d = np.linalg.solve(k_global_bc, global_load_bc)
# print(np.sum(global_load_bc))
t_2_solver = time.perf_counter()
n_dimension = k_global.shape[0]
displm_compelete = np.zeros(n_dimension)
i = 0
j = 0
while i < n_dimension:
if i in bc:
i += 1
else:
displm_compelete[i] = d[j]
i += 1
j += 1
number_lobatto_node = lobatto_pw.shape[0]
number_element_u = len(element_boundaries_u) - 1
number_element_v = len(element_boundaries_v) - 1
number_node_one_row = number_element_u*(number_lobatto_node - 1) + 1
number_node_one_column = number_element_v*(number_lobatto_node - 1) + 1
node_global_a = 1 # u = v = 0 . Four nodes at the tips of the square in u-v parametric space
node_global_b = node_global_a + number_node_one_row - 1
node_global_c = node_global_a + number_element_v*(number_lobatto_node-1)\
*number_node_one_row # u = 0, v = 1
if i_main%2 ==0:
mid_point = int(node_global_a + \
(i_main/2) * number_node_one_row)
else:
mid_point = int(node_global_a + \
((i_main+1)/2) * number_node_one_row)
# print('\nDisplacement ratio zero point: {}'.format(displm_compelete[0]/u_analytic))
# print('\nDisplacement ratio: {}'.format(displm_compelete[5*mid_point - 5]/u_analytic))
# order_displm_array[order_counter] = [i_main, d[5*mid_point - 5]/u_analytic]
print('\nDisplacement ratio: {}'.format(displm_compelete[5*node_global_c - 5]/u_analytic))
order_displm_array[order_counter] = [i_main, displm_compelete[5*node_global_c - 5]/u_analytic]
time_assembling [order_counter] = [i_main, t_2_assembling - t_1_assembling]
time_solver [order_counter] = [i_main, t_2_solver - t_1_solver]
n_dimension_bc = global_load_bc.shape[0] ######################################
dof_displm_array[order_counter] = [n_dimension_bc, \
displm_compelete[5*node_global_c - 5]/u_analytic] ######################################
dof_time_assembling [order_counter] = [n_dimension_bc, t_2_assembling - t_1_assembling] ######################################
dof_time_solver [order_counter] = [n_dimension_bc, t_2_solver - t_1_solver] ######################################
cond_elem [order_counter] = [i_main, np.linalg.cond(k_global_bc)]
order_counter +=1
i_main += 1
np.savetxt(f'curvedbeam_p_ref_displm_elem_uv_{number_element_u}_{number_element_v}.dat', order_displm_array)
np.savetxt(f'curvedbeam_p_ref_asmtime_elem_uv_{number_element_u}_{number_element_v}.dat', time_assembling)
np.savetxt(f'curvedbeam_p_ref_solvertime_elem_uv_{number_element_u}_{number_element_v}.dat', time_solver)
np.savetxt(f'curvedbeam_p_ref_displm_dof_uv_{number_element_u}_{number_element_v}.dat', dof_displm_array)
np.savetxt(f'curvedbeam_p_ref_asmtime_dof_uv_{number_element_u}_{number_element_v}.dat', dof_time_assembling)
np.savetxt(f'curvedbeam_p_ref_solvertime_dof_uv_{number_element_u}_{number_element_v}.dat', dof_time_solver)
np.savetxt(f'curvedbeam_p_ref_cond_elem_uv_{number_element_u}_{number_element_v}.dat', cond_elem)
elif refine_type =="h":
min_order_elem = int(input("\nEnter the minimum order of elements (max order = 30):\n"))
max_order_elem = int(input("Enter the maximum order of elements:\n"))
min_number_elem = int(input("\nEnter the minimum number of elements in u and v direction:\n"))
max_number_elem = int(input("Enter the maximum number of elements in u and v direction:\n"))
print("\nEnter the order of continuity at knots to be used for auto detection of elements boundaries in u direction")
print("The default value is '1'")
c_order_u =int(input())
print("\nEnter the order of continuity at knots to be used for auto detection of elements boundaries in v direction")
print("The default value is '1'")
c_order_v =int(input())
j_main = min_order_elem
while j_main <= max_order_elem:
if j_main==1:
lobatto_pw = lobatto_pw_all[1:3,:]
else:
lobatto_pw_all = lobatto_pw_all # np.delete(lobatto_pw_all,[0, 1, 2], 0)
index = np.argwhere(lobatto_pw_all==j_main)
lobatto_pw = lobatto_pw_all[index[0, 0] + 1:\
index[0, 0] + (j_main+1) + 1, :]
i_main = min_number_elem
elemnum_displm_array = np.zeros((max_number_elem - min_number_elem + 1, 2))
time_assembling = np.zeros((max_number_elem - min_number_elem + 1, 2))
time_solver = np.zeros((max_number_elem - min_number_elem + 1, 2))
dof_displm_array = np.zeros((max_number_elem - min_number_elem + 1, 2)) ######################################
dof_time_assembling = np.zeros((max_number_elem - min_number_elem + 1, 2)) ######################################
dof_time_solver = np.zeros((max_number_elem - min_number_elem + 1, 2)) ######################################
cond_elem =np.zeros((max_number_elem - min_number_elem + 1, 2)) #####################
elemnum_counter = 0
while i_main <= max_number_elem:
print("\n\n\nNumber of elements manually given in u and v: {} Order of elements: {} ".\
format(str(i_main)+'x'+str(i_main), j_main))
print("\nProgram starts to generate mesh according to continuity at knots and manual input of number of elements ...")
u_manual = np.linspace(0, 1, i_main + 1) #np.linspace(a, b, c) divide line bc to c-1 parts or add c points to it.
v_manual = np.linspace(0, 1, i_main + 1)
mesh = gsm.mesh_func(surfs, u_manual, v_manual, c_order_u, c_order_v)
element_boundaries_u = mesh[0]
element_boundaries_v = mesh[1]
number_element_u = len(element_boundaries_u) - 1
number_element_v = len(element_boundaries_v) - 1
bc = gsm.global_boundary_condition(lobatto_pw, bc_h_bott, bc_h_top,\
bc_v_left, bc_v_right, element_boundaries_u,\
element_boundaries_v)
print("\nAssembling global stiffness matrix ...")
t_1_assembling = time.perf_counter()
k_global = gsm.global_stiffness_matrix(surfs, lobatto_pw, \
element_boundaries_u, element_boundaries_v, elastic_modulus, nu)
t_2_assembling = time.perf_counter()
k_global_bc = esm.stiffness_matrix_bc_applied(k_global, bc)
print("\nAssembling global load vector ...")
number_lobatto_node = lobatto_pw.shape[0]
number_element_u = len(element_boundaries_u) - 1
number_element_v = len(element_boundaries_v) - 1
number_node_one_row = number_element_u*(number_lobatto_node - 1) + 1
number_node_one_column = number_element_v*(number_lobatto_node - 1) + 1
node_global_a = 1 # u = v = 0 . Four nodes at the tips of the square in u-v parametric space
node_global_b = node_global_a + number_node_one_row - 1
node_global_c = node_global_a + number_element_v*(number_lobatto_node-1)\
*number_node_one_row # u = 0, v = 1
global_load = glineloadv.global_lineload_vector(surfs, lobatto_pw, \
loaded_side, left_lineload, right_lineload, \
bottom_lineload, top_lineload, element_boundaries_u, \
element_boundaries_v)
global_load_bc = np.delete(global_load, bc, 0)
print("\nLinear solver in action! ...")
t_1_solver = time.perf_counter()
d = np.linalg.solve(k_global_bc, global_load_bc)
t_2_solver = time.perf_counter()
n_dimension = k_global.shape[0]
displm_compelete = np.zeros(n_dimension)
i = 0
j = 0
while i < n_dimension:
if i in bc:
i += 1
else:
displm_compelete[i] = d[j]
i += 1
j += 1
print('\nDisplacement ratio: {}'.format(displm_compelete[5*(node_global_c)-5]/u_analytic))
elemnum_displm_array[elemnum_counter] = [i_main, displm_compelete[5*(node_global_c)-5]/u_analytic] #if number_element_u!=number_element_v then what?????
time_assembling [elemnum_counter] = [i_main, t_2_assembling - t_1_assembling]
time_solver [elemnum_counter] = [i_main, t_2_solver - t_1_solver]
n_dimension_bc = global_load_bc.shape[0] ######################################
dof_displm_array[elemnum_counter] = [n_dimension_bc, \
displm_compelete[5*(node_global_c)-5]/u_analytic] ###################################
dof_time_assembling [elemnum_counter] = [n_dimension_bc, t_2_assembling - t_1_assembling] ######################################
dof_time_solver [elemnum_counter] = [n_dimension_bc, t_2_solver - t_1_solver] #################
# if j_main == max_order_elem:###############
if j_main in [6, 7, 8]:
cond_elem[elemnum_counter] = [i_main, np.linalg.cond(k_global_bc)]
elemnum_counter +=1
i_main += 1
np.savetxt(f'curvedbeam_h_ref_displm_p_{j_main}.dat', elemnum_displm_array)
np.savetxt(f'curvedbeam_h_ref_asmtime_p_{j_main}.dat', time_assembling)
np.savetxt(f'curvedbeam_h_ref_solvertime_p_{j_main}.dat', time_solver)
######################################
np.savetxt(f'curvedbeam_h_ref_displm_dof_p_{j_main}.dat', dof_displm_array)
np.savetxt(f'curvedbeam_h_ref_asmtime_dof_p_{j_main}.dat', dof_time_assembling)
np.savetxt(f'curvedbeam_h_ref_solvertime_dof_p_{j_main}.dat', dof_time_solver)
# if j_main == max_order_elem:###############
if j_main in [6, 7, 8]:##################
np.savetxt(f'curvedbeam_h_ref_cond_elem_p_{j_main}.dat', cond_elem)
j_main += 1
else:
raise ValueError('Refinement can be either "p" or "h"')
# if refine_type == "p" :
# afgnu.plot_displm_solver_assembling_p_ref(number_element_u, number_element_v)
# else:
# afgnu.plot_displm_solver_assembling_h_ref(min_order_elem, max_order_elem)
# # afgnu.plot_displm_time_error_h_ref(min_order_elem, max_order_elem)