### ===================================================================================================================
###  Create NLPO analysis
### ===================================================================================================================
# Copyright ©VIIA 2025
### ===================================================================================================================
###   1. Import modules
### ===================================================================================================================
# General imports
from __future__ import annotations
from typing import TYPE_CHECKING, Optional
# References for functions and classes in the haskoning_structural package
from haskoning_structural.analyses import Analysis
# References for functions and classes in the viiaPackage
if TYPE_CHECKING:
    from viiapackage.viiaStatus import ViiaProject
from viiapackage.analyses.nlpo_helper_functions import viia_center_of_mass_nlpo
from viiapackage.analyses.helper_functions import viia_arclength_control_pivots, viia_find_piles_elements,\
    viia_find_shallow_foundation_elements, viia_find_boundary_interfaces_nodes_elements
from viiapackage.viiaGeneral import viia_find_closest_mesh_node
### ===================================================================================================================
###   2. Function viia_create_analysis_nlpo
### ===================================================================================================================
[docs]def viia_create_analysis_nlpo(
        project: ViiaProject, analysis_nr: str, direction: Optional[str] = None, loadnr: Optional[int] = None,
        ac_pivots: Optional[str] = None, number_eigenmodes: int = 10, modal_pushover_flag: bool = True) -> Analysis:
    """
    Input:
        - project (obj): VIIA project object containing collections of fem objects and project variables.
        - analysis_nr (str): Current analysis referencing the numbers for different analyses. For example 'A4'.
        - direction (str): The directions of which the analysis should be prepared. The default value of the argument is
          None, indicating no directions to be applied. Apply random selection of directions for example:
          'X_neg', 'Y_neg', apply VIIA naming convention for directions.
        - loadnr (int): The associated ID (as integer) of the LoadCombination instance with the name corresponding to
          the NLPO load-cases, if available.
        - ac_pivots (str): Name of method of obtaining arc length control pivots or a sequence of node numbers in string
          format.
        - number_eigenmodes (int): Only used when the modal pushover analysis is set up. This is the number of modes set
          for the eigenvalue analysis executed before the modal pushover analysis, the default value is 10. When the
          dominant mode is after the first 10 modes, this should be adjusted accordingly.
        - modal_pushover_flag (bool): If this argument is set to True, the eigenvalue analysis will executed before
          modal NLPO  analysis, otherwise eigenvalue analysis is not set up for uniform NLPO analysis. The default value
          is True.
    Output:
        - The analysis is created and added to the project.
    """
    # Set analysis name
    if analysis_nr == 'A11':
        name = 'A11 - Niet-lineaire push-over analyse met flexbase - ' + direction
    elif analysis_nr == 'A14':
        name = 'A14 - Niet-lineaire push-over analyse met flexbase met versterkingen - ' + direction
    else:
        raise NotImplementedError("ERROR: Unkown analysis number requested.")
    # Check the supports
    if project.viia_settings.support_type in ['FlexBaseGlobal', 'FlexBaseFinal']:
        pass
    else:
        project.write_log('WARNING: Flex base is not detected for A11 analysis, please check your supports')
    # Initiate list with analysis blocks
    analysis_blocks = list()
    # START This part takes quite long for big model, for now not implemented
    # Get the shapes of horizontal floors, i.e. excluding inclined roof shapes
    # floors = [floor for floor in project.collections.floors
    #           if [round(abs(component), 1) for component in floor.normal_vector()] == [0, 0, 1.0]]
    # floors = [floor for floor in floors if 'FU' not in floor.name and 'N0' not in floor.name]
    #
    # Check if the mesh exists already for the shapes
    # if any(floor.mesh is None for floor in floors):
    #     project.diana_settings.dat_file.get_mesh_from_datfile()
    # END
    # Unless the mesh node is not found or the step before takes too much time
    # Evaluation nodes (SL)
    # Extend evaluation nodes with all floor nodes above ground level
    floor_nodes_dict = project.viia_find_floor_mesh_nodes()
    floor_nodes = [
        node for storey, nodes in floor_nodes_dict.items()
        if 'FU' not in storey and 'N0' not in storey for node in nodes]
    # Add floor nodes to evaluation nodes
    eval_nodes_displa_floors = floor_nodes
    # Find evaluation nodes coordinates from viia_center_of_mass_nlpo()
    eval_coord_dict = viia_center_of_mass_nlpo(project=project)
    eval_coord = [[level['x'], level['y'], level['z']] for level in eval_coord_dict.values()]
    # Get the evaluation nodes according to coordinates
    eval_nodes_displa = [
        viia_find_closest_mesh_node(project=project, target_point=coord, mesh_nodes=floor_nodes)
        for coord in eval_coord]
    # Find the direction object
    direction_obj = None
    if direction == 'X_pos' or direction == 'X_neg':
        direction_obj = project.get_direction(name='X')
    elif direction == 'Y_pos' or direction == 'Y_neg':
        direction_obj = project.get_direction(name='Y')
    # Define the arclength control points based om default procedure or provided ac_pivots
    if ac_pivots:
        pivots_by_method = viia_arclength_control_pivots(project, ac_pivots)
        if pivots_by_method[0] is None:
            project.write_log("WARNING: Attempting to interpret the input for arclength control pivots as "
                              "predefined node numbers.", True)
            pivot_nodes = ac_pivots
        else:
            pivot_nodes = pivots_by_method
    else:
        project.write_log("WARNING: SDF pushover analysis requires definition of arclength control pivots. Default "
                          "method is being called.", True)
        pivot_nodes = viia_arclength_control_pivots(project)
    # Uniform or modal pushover
    # Adding specific eigen mode for different analysis direction
    x_mode = project.project_specific['modes'][0]
    y_mode = project.project_specific['modes'][1]
    dominant_mode = None
    if direction == "X_pos" or direction == 'X_neg':
        dominant_mode = x_mode
        # project.diana_settings.analysis_settings['defaults_viia']['EIGEN']['output'][
        #     'SXXX_OUTPUT_EIGENVALUE'].update({'modes': str(dominant_mode)})
    if direction == 'Y_pos' or direction == 'Y_neg':
        dominant_mode = y_mode
        # project.diana_settings.analysis_settings['defaults_viia']['EIGEN']['output'][
        #     'SXXX_OUTPUT_EIGENVALUE'].update({'modes': str(dominant_mode)})
    if modal_pushover_flag:
        # Create the eigenvalue analysis as part of the NLPO
        execute_block = project.create_specified_execute_eigenvalue(maximum_number_iterations=30)
        output_block = list()
        output_block.append(project.create_output_block_eigenvalue(
            name='OUTPUT_EIGENVALUE',
            output_filename=project.name + '_v' + str(project.version) + '_OUTPUT_EIGENVALUE',
            modes=[dominant_mode],
            manual_nodes=eval_nodes_displa_floors,
            output_device='tabulated'
        ))
        # Create the analysis blocks for the eigenvalue analysis (in this case one block only)
        analysis_blocks.append(project.create_eigenvalue_analysis_block(
            name='Structural eigenvalue', modes=number_eigenmodes, execute=execute_block, output_blocks=output_block))
    # Create the execute blocks for the NLPO
    analysis_block_name = 'Nonlinear Pushover'
    physical_nonlinear_options = project.create_physical_nonlinear_options()
    arclength_settings = project.create_arc_length_control_settings(
        unloading_determination='negative pivot',
        control_sets_node_numbers=pivot_nodes,
        control_set_direction=direction_obj)
    # Obtaining the nodes for the center of mass on a storey level
    project.viia_get_drift_limits(eval_coord)
    stop_criteria = []
    for i in range(len(project.project_specific['drift_limits'])):
        if 'x' in direction.lower():
            displacement_stop_criteria = project.create_output_item('displacements', component='x')[0]
        elif 'y' in direction.lower():
            displacement_stop_criteria = project.create_output_item('displacements', component='y')[0]
        else:
            raise ValueError(f"ERROR: stop criteria component could not be determined for direction {direction}.")
        if i == 0:
            stop_criteria.append(project.create_stop_criteria(
                selected_nodes=[eval_nodes_displa[i]], selected_elements='all',
                displacement_max_value=project.project_specific['drift_limits']['N'+str(i)],
                displacement_min_value=-1*project.project_specific['drift_limits']['N'+str(i)],
                displacement_stop_criteria=displacement_stop_criteria))
        else:
            stop_criteria.append(
                project.create_stop_criteria(
                    selected_nodes=[eval_nodes_displa[i]],
                    displacement_max_value=project.project_specific['drift_limits']['N' + str(i)],
                    displacement_min_value=-1*project.project_specific['drift_limits']['N'+str(i)],
                    displacement_stop_criteria=displacement_stop_criteria,
                    base_node=eval_nodes_displa[i-1]))
    execute_blocks = [
        project.create_specified_execute_nonlinear(
            name='SW',
            physical_nonlinear_options=None,
            satisfy_all_criteria=True,
            line_search=False,
            steps=project.create_load_steps(load_combination=project.viia_get('load_combinations', name='DL')),
            maximum_number_iterations=30,
            convergence_displacement=True,
            convergence_displacement_tolerance=0.001,
            convergence_displacement_abort=1.0E+17,
            convergence_force=True,
            convergence_force_tolerance=0.001,
            convergence_force_abort=1.0E+17,
            iteration_method='newton-raphson',
            continue_iteration=True),
        project.create_specified_execute_nonlinear(
            physical_nonlinear_options=physical_nonlinear_options,
            line_search=True,
            name='FIRST PUSH',
            steps=project.create_load_steps(
                steps=[[0.1, 1]], load_combination=project.get_load_combination(_id=loadnr)),
            maximum_number_iterations=30,
            convergence_displacement=True,
            convergence_displacement_tolerance=0.01,
            convergence_displacement_abort=1.0E+4,
            iteration_method='newton-raphson',
            continue_iteration=True),
        project.create_specified_execute_nonlinear(
            physical_nonlinear_options=None,
            line_search=True,
            continue_iteration=True,
            name='NLPO',
            steps=project.create_load_steps(
                steps=[[0.1, 500]], load_combination=project.viia_get('load_combinations', id=loadnr),
                arclength_control=arclength_settings),
            maximum_number_iterations=30,
            convergence_displacement=True,
            convergence_displacement_tolerance=0.01,
            convergence_displacement_abort=1.0E+17,
            iteration_method='newton-raphson',
            stop_criterion=stop_criteria)]
    # Flexbase analysis (SL) - This part is equal for all pushovers
    # 1. Collect nodes and elements
    # 2. Set the output block into traction (for shallow foundation) / nodal force (pile foundation)
    output_blocks = list()
    output_name = 'OUTPUT_NLPO_SL_NODES'
    output_filename = project.name + '_v' + str(project.version) + '_' + output_name
    location_option = project.create_element_output_item_options(location='integration points')
    output_blocks.append(project.create_output_block_nonlinear(
        output_device=None,
        output_filename=output_filename,
        output_items=[
            project.create_output_item(output='displacements'),
            project.create_output_item(
                output='strain', output_type='crack width', theoretical_formulation='green-lagrange',
                operation='principal', options=location_option),
            project.create_output_item(
                output='stress', theoretical_formulation='cauchy', options=location_option),
            project.create_output_item(
                output='strain', output_type='total', theoretical_formulation='relative displacement',
                operation='local', options=location_option),
            project.create_output_item(
                output='stress', theoretical_formulation='cauchy', operation='principal', options=location_option),
            project.create_output_item(output='reaction forces'),
            project.create_output_item(output='residual forces')],
        name=output_name,
        steps='All'))
    output_name = 'OUTPUT_NLPO_SL_INTPNT'
    output_filename = project.name + '_v' + str(project.version) + '_' + output_name
    output_blocks.append(project.create_output_block_nonlinear(
        output_device=None,
        output_filename=output_filename,
        output_items=[
            project.create_output_item(
                output='stress',
                output_type='total',
                theoretical_formulation='traction',
                operation='local',
                options=location_option)],
        name=output_name,
        steps='All'))
    # Setting TB file (evaluations nodes disp v.s. base shear)
    # Evaluation point
    output_name = 'OUTPUT_NLPO_DISPLA_TB'
    output_filename = project.name + '_v' + str(project.version) + '_' + output_name
    output_blocks.append(project.create_output_block_nonlinear(
        output_device='tabulated',
        output_filename=output_filename,
        output_items=[project.create_output_item('displacements')],
        name=output_name,
        manual_nodes=eval_nodes_displa_floors,
        steps='All'))
    # Create output blocks for NLPO if it is only pile or shallow foundations
    # Only pile foundation
    mixed_foundation = False
    location_option = project.create_element_output_item_options(location='integration points')
    if project.collections.piles != [] and project.collections.boundary_interfaces == []:
        # Collection nodes and elements of the spring on the pile
        eval_nodes_base_shear, eval_elements_base_shear = project.viia_find_top_huanbeam_nodes_elements()
        output_item_base_shear = [
            project.create_output_item(
                output='nodal element force',
                output_type='total',
                theoretical_formulation='translation',
                operation='global')]
    # Only shallow foundation
    elif project.collections.boundary_interfaces != [] and project.collections.piles == []:
        # Collection of node and element on the boundary interfaces for shallow foundation
        eval_nodes_base_shear, eval_elements_base_shear = viia_find_boundary_interfaces_nodes_elements(project=project)
        output_item_base_shear = [
            project.create_output_item('reaction forces'),
            project.create_output_item(
                output='stress',
                theoretical_formulation='traction',
                operation='local',
                options=location_option)]
    # Mixed
    # First collect data of piles, then foundations
    else:
        mixed_foundation = True
        eval_elements_base_shear, eval_nodes_base_shear, _, _, _ = viia_find_piles_elements(project=project)
        output_item_base_shear = [
            project.create_output_item(
                output='nodal element force',
                output_type='total',
                theoretical_formulation='translation',
                operation='global')]
    output_name = 'OUTPUT_NLPO_BASE_SHEAR_TB'
    output_filename = project.name + '_v' + str(project.version) + '_' + output_name
    output_blocks.append(project.create_output_block_nonlinear(
        output_device='tabulated',
        output_filename=output_filename,
        output_items=output_item_base_shear,
        name=output_name,
        manual_nodes=eval_nodes_base_shear,
        manual_elements=eval_elements_base_shear,
        steps='All'))
    # Add extra block for mixed foundation
    if mixed_foundation:
        output_name = 'OUTPUT_NLPO_BASE_SHEAR_TB_fstrips'
        output_filename = project.name + '_v' + str(project.version) + '_' + output_name
        output_item_base_shear = [
            project.create_output_item(
                output='reaction forces'),
            project.create_output_item(
                output='stress',
                theoretical_formulation='traction',
                operation='local',
                options=location_option)]
        shallow_foundation_elements, shallow_foundation_nodes, _, _ = viia_find_shallow_foundation_elements(project)
        output_blocks.append(project.create_output_block_nonlinear(
            output_device='tabulated',
            output_filename=output_filename,
            output_items=output_item_base_shear,
            name=output_name,
            manual_nodes=shallow_foundation_nodes,
            manual_elements=shallow_foundation_elements,
            steps='All'))
    # Set up the evaluation model
    evaluation_model = project.create_specified_evaluate_model_settings()
    # Set up the NLPO settings
    physical_nonlinearity = project.create_specified_physical_nonlinearity()
    geometrical_nonlinearity = project.create_specified_geometrical_nonlinearity(
        geometrical_nonlinearity_type='Total Lagrange')
    nlpo_settings = project.create_specified_nonlinear_settings(
        physical_nonlinearity=physical_nonlinearity,
        geometrical_nonlinearity=geometrical_nonlinearity,
        transient_effects=None,
        linearisation=False,
        recomputation=False)
    analysis_blocks.append(project.create_nonlinear_analysis_block(
        execute_blocks=execute_blocks,
        output_blocks=output_blocks,
        name=analysis_block_name,
        evaluate_model=evaluation_model,
        nonlinear_settings=nlpo_settings))
    # Creation of the analysis object
    return project.create_analysis(name=name, calculation_blocks=analysis_blocks) 
### ===================================================================================================================
###   3. End of script
### ===================================================================================================================