### ===================================================================================================================
### Function to calculate the base shear
### ===================================================================================================================
# Copyright ©VIIA 2024
### ===================================================================================================================
### 1. Import modules
### ===================================================================================================================
# General imports
from __future__ import annotations
from typing import TYPE_CHECKING, List, Union, Optional
from pathlib import Path
import json
# References for functions and classes in the rhdhv_fem package
from rhdhv_fem.fem_tools import fem_make_list_of_text, fem_read_file
from rhdhv_fem.analyses.analysis import Analysis
from rhdhv_fem.fem_config import Config
from rhdhv_fem.fem_math import fem_greater, fem_smaller
# References for functions and classes in the viiaPackage
if TYPE_CHECKING:
from viiapackage.viiaStatus import ViiaProject
# Import module matplotlib
import matplotlib
# Switch to a non-interactive backend to properly import matplotlib.pyplot
matplotlib.use(Config.MPL_NONINTERACTIVE(notify=False))
import matplotlib.pyplot as plt
### ===================================================================================================================
### 2. Function viia_base_shear
### ===================================================================================================================
[docs]def viia_base_shear(
project: ViiaProject, analysis: Analysis, base_shear_tb_file: List[Union[Path, str]] = None,
geo_postprocessed_data_json_file_abaqus: Union[Path, str] = None):
"""
This function calculates the base shear from a DIANA analysis or ABAQUS geo post-processing for an object with
strip foundation, or pile foundation or mixed foundations.
Input:
- project (obj): Project object containing collections and of fem objects and project variables.
- base_shear_tb_file(list): List of output .tb files for calculating the base shear from DIANA.
- geo_postprocessed_data_json_file_abaqus(string or Path): Path or string of path of output .json
files for geo postprocessing from ABAQUS.
- analysis (obj): Object reference of the analysis of which the base shear should be calculated.
Output:
- A figure with the base shear in time-history is created and saved in the current analysis folder.
- A json file with calculated base shear in time-history is created and saved in the current analysis folder.
"""
def plot_base_shear(
title: str, file_name: str, time_data: list, x_data: list, y_data: list, z_data: list, max_x: Optional[str],
max_y: Optional[str], max_z: Optional[str] = None, min_z: Optional[str] = None):
""" Function to plot the base shear figures."""
# Create graph with base shear results
plt.figure(figsize=(15, 5))
# Apply VIIA graph style
style_file = Path(project.viia_settings.project_specific_package_location) / 'viiaGraph.mplstyle'
plt.style.use(style_file.as_posix())
# Plot the values
fig1, ax1 = plt.subplots(figsize=(15, 5))
ax1.plot(time_data, x_data, label='Total base shear in x-direction')
ax1.plot(time_data, y_data, label='Total base shear in y-direction')
# Adjust the axes limits
y_min = min(min(x_data) * 1.05, min(y_data) * 1.05)
y_max = max(max(x_data) * 1.05, max(y_data) * 1.05)
plt.ylim(y_min, y_max)
# Apply formatting
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('Base shear [kN]')
ax1.xaxis.set_ticks_position('bottom')
if max_x and max_y:
ax1.set_title(f"{title} (Maximum base shear: X {max_x} kN, Y {max_y} kN)")
else:
ax1.set_title(title)
ax1.legend(loc='upper right')
# Save the horizontal base shear graph
fig1.savefig(project.current_analysis_folder.as_posix() + f'/{file_name}' + '_Horizontal.png')
plt.close()
# Create graph with vertical base shear results
plt.figure(figsize=(15, 5))
# Apply VIIA graph style
style_file = Path(project.viia_settings.project_specific_package_location) / 'viiaGraph.mplstyle'
plt.style.use(style_file.as_posix())
# Plot the values
fig2, ax2 = plt.subplots(figsize=(15, 5))
ax2.plot(time_data, z_data, label='Total vertical load')
# Adjust the axes limits
y_min = min(z_data) * 1.05
y_max = max(z_data) * 1.05
plt.ylim(y_min, y_max)
# Apply formatting
ax2.set_xlabel('Time [s]')
ax2.set_ylabel('Base shear [kN]')
ax2.xaxis.set_ticks_position('bottom')
if max_z and min_z:
ax2.set_title(f"{title} (Maximum vertical force: {max_z} kN, Minimum vertical force: {min_z} kN)")
else:
ax2.set_title(title)
ax2.legend(loc='upper right')
# Save the vertical base shear graph
fig2.savefig(project.current_analysis_folder.as_posix() + f'/{file_name}' + '_Vertical.png')
plt.close()
# Check which files are provided
if base_shear_tb_file is None and geo_postprocessed_data_json_file_abaqus is None:
raise ValueError('ERROR: No files have been selected. Please check the input.')
elif base_shear_tb_file is not None and geo_postprocessed_data_json_file_abaqus is not None:
raise ValueError('ERROR: Only one type of files should be selected. Please check the input.')
# Get analysis number
analyse_nr = analysis.name.split(' - ')[0]
# Check if shallow foundation is present and if piles are present
fos = project.shallow_foundation_present()
piles = project.pile_foundation_present()
if not piles and not fos:
raise ValueError(
"ERROR: No shallow or pile foundation present in the model, please check. Maybe you are using the model "
"json not the analysis json.")
# If tb-file is provided:
if base_shear_tb_file is not None:
if not isinstance(base_shear_tb_file, list):
base_shear_tb_file = [base_shear_tb_file]
for i, file in enumerate(base_shear_tb_file):
if isinstance(file, str):
base_shear_tb_file[i] = Path(file)
elements = {}
# Procedure for shallow foundations
if fos:
# Read file
fos_file = None
for file in base_shear_tb_file:
if file is not None:
if '5A.tb' in file.as_posix():
fos_file = file
if not fos_file:
project.write_log("ERROR: Input file not recognized for base shear function, check your input.")
return
lines = fem_read_file(fos_file)
if not lines:
return
if project.viia_settings.support_type == 'FixedBase':
time_bool = False
current_time = 0
base_shear_x_fos = []
base_shear_y_fos = []
base_shear_z_fos = []
time_fos = []
for i, line in enumerate(lines):
line = fem_make_list_of_text(line)
if 'Time' in line:
if not time_bool:
time_bool = True
else:
time_fos.append(current_time)
base_shear_x_fos.append(base_shear_x_fos_timestep)
base_shear_y_fos.append(base_shear_y_fos_timestep)
base_shear_z_fos.append(base_shear_z_fos_timestep)
base_shear_x_fos_timestep = 0
base_shear_y_fos_timestep = 0
base_shear_z_fos_timestep = 0
current_time = float(line[1])
if len(line) == 4 and time_bool:
try:
int(line[0])
base_shear_x_fos_timestep += float(line[1])
base_shear_y_fos_timestep += float(line[2])
base_shear_z_fos_timestep += float(line[3])
except ValueError:
pass
else:
mesh_elem_area = {}
for interface in project.collections.boundary_interfaces:
for elem in interface.mesh.mesh_elements:
mesh_elem_area[elem.id] = {'area': elem.get_area()}
time_bool = False
current_time = 0
current_element = None
base_shear_x_fos = []
base_shear_y_fos = []
base_shear_z_fos = []
time_fos = []
element_sx = 0
element_sy = 0
element_nz = 0
nr_int_points = 1
base_shear_x_fos_timestep = 0
base_shear_y_fos_timestep = 0
base_shear_z_fos_timestep = 0
for i, line in enumerate(lines):
line = fem_make_list_of_text(line)
if 'Time' in line:
if not time_bool:
time_bool = True
else:
time_fos.append(current_time)
base_shear_x_fos.append(base_shear_x_fos_timestep)
base_shear_y_fos.append(base_shear_y_fos_timestep)
base_shear_z_fos.append(base_shear_z_fos_timestep)
base_shear_x_fos_timestep = 0
base_shear_y_fos_timestep = 0
base_shear_z_fos_timestep = 0
current_time = float(line[1])
elif len(line) == 5 and time_bool:
try:
new_element = int(line[0])
if elements == {} or new_element not in elements:
elements[new_element] = {}
if current_element:
base_shear_x_fos_timestep += \
float(element_sx) / nr_int_points * mesh_elem_area[current_element]['area']
if elements[new_element] == {}:
elements[new_element]['base_shear_x'] = []
elements[new_element]['base_shear_x'].append(
float(element_sx) / nr_int_points * mesh_elem_area[current_element]['area'])
base_shear_y_fos_timestep += \
float(element_sy) / nr_int_points * mesh_elem_area[current_element]['area']
if 'base_shear_y' not in elements[new_element]:
elements[new_element]['base_shear_y'] = []
elements[new_element]['base_shear_y'].append(
float(element_sy) / nr_int_points * mesh_elem_area[current_element]['area'])
base_shear_z_fos_timestep += \
float(element_nz) / nr_int_points * mesh_elem_area[current_element]['area']
if 'base_shear_z' not in elements[new_element]:
elements[new_element]['base_shear_z'] = []
elements[new_element]['base_shear_z'].append(
float(element_nz) / nr_int_points * mesh_elem_area[current_element]['area'])
current_element = new_element
element_sx = float(line[2])
element_sy = float(line[3])
element_nz = float(line[4])
nr_int_points = 1
except ValueError:
pass
elif len(line) == 4 and time_bool:
try:
element_sx += float(line[1])
element_sy += float(line[2])
element_nz += float(line[3])
nr_int_points += 1
except ValueError:
pass
# If pile foundation present
if piles:
# Read file
piles_file = None
for file in base_shear_tb_file:
if file is not None:
if '5B.tb' in file.as_posix():
piles_file = file
if not piles_file:
project.write_log("ERROR: Input file not recognised for base shear function, check your input.")
return
lines = fem_read_file(piles_file)
if not lines:
return
time_bool = False
node_force = False
current_time = 0
base_shear_x_pile = []
base_shear_y_pile = []
base_shear_z_pile = []
time_pile = []
for i, line in enumerate(lines):
line = fem_make_list_of_text(line)
if 'Time' in line:
next_line = fem_make_list_of_text(lines[i + 1])
if 'DISPLA' in next_line:
node_force = False
continue
elif 'NODFOR' in next_line:
node_force = True
if not time_bool:
time_bool = True
else:
time_pile.append(current_time)
base_shear_x_pile.append(base_shear_x_pile_timestep)
base_shear_y_pile.append(base_shear_y_pile_timestep)
base_shear_z_pile.append(base_shear_z_pile_timestep)
base_shear_x_pile_timestep = 0
base_shear_y_pile_timestep = 0
base_shear_z_pile_timestep = 0
current_time = float(line[1])
if len(line) == 4 and time_bool and node_force:
try:
base_shear_x_pile_timestep += float(line[1])
base_shear_y_pile_timestep += float(line[2])
base_shear_z_pile_timestep += float(line[3])
current_element = int(line[0])
if elements == {} or current_element not in elements:
elements[current_element] = {}
if 'base_shear_x' not in elements[current_element]:
elements[current_element]['base_shear_x'] = []
elements[current_element]['base_shear_x'].append(float(line[1]))
if 'base_shear_y' not in elements[current_element]:
elements[current_element]['base_shear_y'] = []
elements[current_element]['base_shear_y'].append(float(line[2]))
if 'base_shear_z' not in elements[current_element]:
elements[current_element]['base_shear_z'] = []
elements[current_element]['base_shear_z'].append(float(line[3]))
except ValueError:
pass
if fos and not piles:
base_shear_x = base_shear_x_fos
base_shear_y = base_shear_y_fos
base_shear_z = base_shear_z_fos
time = time_fos
elif piles and not fos:
base_shear_x = base_shear_x_pile
base_shear_y = base_shear_y_pile
base_shear_z = base_shear_z_pile
time = time_pile
else:
if len(time_fos) != len(time_pile):
project.write_log(
"ERROR: Number of timesteps in FOS output does not correspond to " +
"the number of timesteps in pile output. Check files.", True)
return
time = time_fos
base_shear_x = []
base_shear_y = []
base_shear_z = []
for i in range(len(time)):
base_shear_x.append(base_shear_x_fos[i] + base_shear_x_pile[i])
base_shear_y.append(base_shear_y_fos[i] + base_shear_x_pile[i])
base_shear_z.append(base_shear_z_fos[i] + base_shear_x_pile[i])
# Get signal name
if base_shear_tb_file[0]:
signal_name = base_shear_tb_file[0].parts[-2]
elif base_shear_tb_file[1]:
signal_name = base_shear_tb_file[1].parts[-2]
if analyse_nr not in project.results:
project.results[analyse_nr] = {}
# Prepare dump file
name = 'base_shear.json'
if base_shear_tb_file[0]:
dumpfile = base_shear_tb_file[0].parent / name
else:
dumpfile = base_shear_tb_file[1].parent / name
# if geo_postprocessed_data.json is provided
if geo_postprocessed_data_json_file_abaqus is not None:
file = geo_postprocessed_data_json_file_abaqus
if isinstance(file, str):
geo_postprocessed_data_json_file_abaqus = Path(file)
# Read file
file = geo_postprocessed_data_json_file_abaqus
with open(file) as json_file:
s_all = json.load(json_file)
for item in s_all['section_force_requests']['items']:
if "force" not in item.keys():
project.write_log("ERROR: Base shear forces are not recognized in the input file, check your input.")
return
time = list(map(
float, s_all['section_force_requests']['items'][0]['force']['1']['Structural nonlinear 2'].keys()))
# Get time history from shear force time history
base_shear_x = []
base_shear_y = []
base_shear_z = []
for time_step in time:
force_x = 0
force_y = 0
force_z = 0
for item in s_all['section_force_requests']['items']:
force_x += item['force']['1']['Structural nonlinear 2'][str(time_step)]
force_y += item['force']['2']['Structural nonlinear 2'][str(time_step)]
force_z += item['force']['3']['Structural nonlinear 2'][str(time_step)]
base_shear_x.append(force_x)
base_shear_y.append(force_y)
base_shear_z.append(-force_z)
elements = {}
# Get signal name
if '/' in geo_postprocessed_data_json_file_abaqus.as_posix():
signal_name = geo_postprocessed_data_json_file_abaqus.name.split('_')[0].split('-')[-1]
else:
signal_name = geo_postprocessed_data_json_file_abaqus.name.split('_')[0].split('-')[-1]
if analyse_nr not in project.results:
project.results[analyse_nr] = {}
# Prepare dump file
name = f'{signal_name}_Base_Shear_Input_NLTH.json'
dumpfile = geo_postprocessed_data_json_file_abaqus.parent / name
# Logging min and maximum base shear in all direction.
max_base_shear_x = max(base_shear_x)
min_base_shear_x = min(base_shear_x)
max_base_shear_y = max(base_shear_y)
min_base_shear_y = min(base_shear_y)
max_base_shear_z = max(base_shear_z)
min_base_shear_z = min(base_shear_z)
# Logging timestep of min and maximum base shear in all direction
crit_timesteps = dict()
crit_timesteps['max_base_shear_x_time'] = base_shear_x.index(max_base_shear_x) + 1
crit_timesteps['max_base_shear_y_time'] = base_shear_y.index(max_base_shear_y) + 1
crit_timesteps['max_base_shear_z_time'] = base_shear_z.index(max_base_shear_z) + 1
crit_timesteps['min_base_shear_z_time'] = base_shear_z.index(min_base_shear_z) + 1
# Get absolute value of maximum base shear in horizontal direction
max_base_shear_x = max(abs(min_base_shear_x), max_base_shear_x)
max_base_shear_y = max(abs(min_base_shear_y), max_base_shear_y)
# Add to results dictionary
project.results[analyse_nr]['base_shear'] = {}
project.results[analyse_nr]['base_shear'][signal_name] = {
'base_shear_x': base_shear_x,
'base_shear_y': base_shear_y,
'vertical_force_z': base_shear_z,
'time': time,
'max_base_shear_x': max_base_shear_x,
'max_base_shear_y': max_base_shear_y,
'max_normal_force_z': max_base_shear_z,
'min_normal_force_z': min_base_shear_z,
'force_element_time': elements,
'crit_timesteps': crit_timesteps}
# Logging min and maximum base shear in all direction.
max_base_shear_x = str(round(abs(max_base_shear_x) / 1000, 0))
max_base_shear_y = str(round(abs(max_base_shear_y) / 1000, 0))
max_normalforce_z = str(round(max_base_shear_z / 1000, 0))
min_normalforce_z = str(round(min_base_shear_z / 1000, 0))
project.write_log(
f"The maximum base shear in x-direction is {max_base_shear_x} kN at timestep "
f"{crit_timesteps['max_base_shear_x_time']}.")
project.write_log(
f"The maximum base shear in y-direction is {max_base_shear_y} kN at timestep "
f"{crit_timesteps['max_base_shear_y_time']}.")
project.write_log(
f"The maximum total normal force in z-direction is {max_normalforce_z} kN at timestep "
f"{crit_timesteps['max_base_shear_z_time']}.")
project.write_log(
f"The minimum total normal force in z-direction is {min_normalforce_z} kN at timestep "
f"{crit_timesteps['min_base_shear_z_time']}.")
# Revert values from [N] to [kN] for visualisation
for i in range(len(base_shear_x)):
base_shear_x[i] = base_shear_x[i] / 1000
for i in range(len(base_shear_y)):
base_shear_y[i] = base_shear_y[i] / 1000
for i in range(len(base_shear_z)):
base_shear_z[i] = base_shear_z[i] / 1000
# If pile and shallow foundation are combined
if piles and fos:
for i in range(len(base_shear_x_pile)):
base_shear_x_pile[i] = base_shear_x_pile[i] / 1000
base_shear_x_fos[i] = base_shear_x_fos[i] / 1000
for i in range(len(base_shear_y_pile)):
base_shear_y_pile[i] = base_shear_y_pile[i] / 1000
base_shear_y_fos[i] = base_shear_y_fos[i] / 1000
for i in range(len(base_shear_z_pile)):
base_shear_z_pile[i] = base_shear_z_pile[i] / 1000
base_shear_z_fos[i] = base_shear_z_fos[i] / 1000
# Create picture with the total base shear
plot_base_shear(
title=f"Base shear analysis {analysis.name.split(' -')[0]} with {signal_name}", file_name='BaseShear',
time_data=time, x_data=base_shear_x, y_data=base_shear_y, z_data=base_shear_z, max_x=max_base_shear_x,
max_y=max_base_shear_y, max_z=max_normalforce_z, min_z=min_normalforce_z)
if fos and piles:
# fos base shear
plot_base_shear(
title=f"Base shear analysis {analysis.name.split(' -')[0]} with {signal_name}, shallow foundation",
file_name='BaseShear_fos', time_data=time, x_data=base_shear_x_fos, y_data=base_shear_y_fos,
z_data=base_shear_z_fos, max_x=None, max_y=None)
# pile base shear
plot_base_shear(
title=f"Base shear analysis {analysis.name.split(' -')[0]} with {signal_name}, pile foundation",
file_name='BaseShear_piles', time_data=time, x_data=base_shear_x_pile, y_data=base_shear_y_pile,
z_data=base_shear_z_pile, max_x=None, max_y=None)
# Collect the data for the base signal graph
signal_t = []
signal_x = []
signal_y = []
signal_z = []
for timeseries_combination in project.collections.timeseries_combinations:
if timeseries_combination.name.split('_')[-1] == signal_name:
signal_x = timeseries_combination.x_direction.value_series
signal_y = timeseries_combination.y_direction.value_series
signal_z = timeseries_combination.z_direction.value_series
signal_t = timeseries_combination.x_direction.time_series
# Create the base signal graph
plt.figure(figsize=(15, 5))
viia_style_file = Path(project.viia_settings.project_specific_package_location) / 'viiaGraph.mplstyle'
plt.style.use(viia_style_file.as_posix())
ax1 = plt.subplot2grid((1, 1), (0, 0))
ax1.plot(signal_t, [sx / project.gravitational_acceleration for sx in signal_x], label='Base signal in x-direction')
ax1.plot(signal_t, [sy / project.gravitational_acceleration for sy in signal_y], label='Base signal in y-direction')
ax1.plot(signal_t, [sz / project.gravitational_acceleration for sz in signal_z], label='Base signal in z-direction')
ax1.spines['bottom'].set_position('zero')
plt.xlabel('Time [s]')
plt.ylabel('Acceleration [g]')
plt.title('Base signal ' + signal_name)
plt.legend(loc='upper right')
plt.savefig(project.current_analysis_folder.as_posix() + f"/BaseSignal_{signal_name}.png")
plt.close()
with open(dumpfile, 'w') as fd:
json.dump(project.results[analyse_nr], fd, indent=2, sort_keys=True)
return
### ===================================================================================================================
### 3. End of script
### ===================================================================================================================