# -*- coding: utf-8 -*-
"""
Created on Fri Aug 25 13:08:16 2020

@author: haolinl
"""

import copy
import os
import time

import numpy as np
import random
import scipy.io # For extracting data from .mat file


class inputFileGenerator(object):
    """
    Generate input file for Abaqus. 

    Unit system: 
        Length: m
        Force: N
        Pressure: Pa
    """

    def __init__(self, data_file_name, write_path, material_type, fix_indices_list, node_variable_name, elem_variable_name, user_prescribed_force_field=[]):
        """
        Initialize parameters. 

        Parameters: 
        ----------
        data_file_name: String. 
            The file path of information of node, element, etc. 
        write_path: String. 
            The path to write the inp file. 
        material_type: String. 
            The type of material. 
            Used to indicate whether to consider material nonlinearity.
        fix_indices_list: List of ints. 
            The node indices to be fixed. 
        node_variable_name: String. 
            The variable name of the nodes matrix in the data file. 
        elem_variable_name: String. 
            The variable name of the elements matrix in the data file. 
        user_prescribed_force_field (optional): List of floats.
            The user-prescribed vector of laplacian force field. 
            Size: nSurfI x 3. 
            Default: []. 
        """
        
        # Data & Variables. 
        self.data_file_name = data_file_name
        self.data_mat = scipy.io.loadmat(self.data_file_name)
        self._surface_mat = self.data_mat["FaceI"]
        self._surface_nodes = self.data_mat["idxSurfI"]
        self._surface_nodes_num = self.data_mat["nSurfI"][0,0]
        self._outer_surface_regionNum = 22 # Int. The region number of outer surface. 
        self._outer_surface_nodes_list = self._extractOuterSurfaceNodes(self.data_mat["faces"], self._outer_surface_regionNum) # List of sorted ints. The indices of outer surface nodes. Indexed from 1. 
        self._outer_surface_nodes_num = len(self._outer_surface_nodes_list)
        self._triangle_nodes_list = []
        self._coupled_list = []
        self._node_variable_name = node_variable_name
        self._elem_variable_name = elem_variable_name
        self._inputFile_lines_total = []
        self.writePath = write_path

        self._modulus = 1e7 # Young's modulus. Unit: Pa. Default: 1e7. 
        self._poisson_ratio = 0.48 # Poisson's ratio. Linear elastic default: 0.3; neo-Hookean default: 0.48.  

        self._isCoupleOn = False # Boolean. True: use coupling constraint; False: do not use coupling constraint. Must not turn on if applying Laplacian smoothing.  
        self._coupling_type = "Kinematic" # String. "Kinematic" / "Distributing". 
        self._coupling_neighbor_layers = 1 # How deep does the neighborhood searching go. Default: 1. 

        self._isLaplacianSmoothingOn = True # Boolean. True: use laplacian smoothing. False: do not use laplacian smoothing.
        self._laplacian_variable_name = "laplacianMatrixI3" 
        self._massMatrix_variable_name = "massMatrixI3" 
        self._laplacian_iter_num = 20 # Default: 3. 
        self._smoothing_rate = 0.1 # Default: 0.1 (Previous: 1e-4). 

        self.loads_num = 3 # For initial testing.
        self._load_sampling_style = "gaussian" # String. Indicating the type of random sampling for force components. "uniform" / "gaussian". 
        self._load_scale = (0.0, 10.0) # Absolute range of the force for uniform sampling. Case and BC specific. (min, max). Unit: N.
        self._gaussian_params = (4.0, 0.8) # Mean and deviation of the force for Gaussian sampling. Case and BC specific. (mean, deviation). Unit: N.
        self._load_params_tuple = None
        self._initial_force_component_vector = [] # List of floats. Default: []. Example: [5., 5., 5.]. 

        self.autoIncrementNum = 5000 # Int. The maximum increment number of the AutoSolver. 
        self.initIncrem = 0.001 # Float. The initial length of the increment (for fixed-step, this is also the length per increm). 
        self.minIncrem = 1e-20 # Float. The minimum increment length for the AutoSolver (ueless for the StaticSolver). 
        self.maxIncrem = 1.0 # Float. The maximum increment length for the AutoSolver (useless for the StaticSovler). 
        self.totalTime = 1.0 # Float. The total time for one simulation step. 
        self.frameNum = 1 # Int. The number of frames intending to extract from the nodal file.

        # ================== Load sampling variables ================== #

        if self._isCoupleOn: self._couple_region_num = self.loads_num
        else: self._couple_region_num = 0

        if self._load_sampling_style == "gaussian": self._load_params_tuple = self._gaussian_params
        elif self._load_sampling_style == "uniform": self._load_params_tuple = self._load_scale
        else: 
            self._load_sampling_style = "uniform"
            self._load_params_tuple = self._load_scale

        # ============================================================= #

        # Header. 
        self._header = ["*Heading"]

        # Part definition. 
        self._part_name = "part-1"
        self._material_name = "tissue"
        self._part_initial = ["*Part, name={}".format(self._part_name)] # Total list of Part definition. 
        self._node = ["*Node"]
        self._elem = ["*Element, type=C3D10"] # Nonlinear tetrahedron. http://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node33.html#tennode. 
        self._nset_all = []
        self._elset_all = []
        self._section = ["*Solid Section, elset=allElems, material={}".format(self._material_name),
                         ","]
        self._part_end = ["*End Part"]
        self._new_node_list = []
        self._new_node_dict = {}
        self._node_num = None
        self._orig_node_num = None
        self._elem_num = None
        self._part = self.generatePart()

        # Load settings.
        self._loads_nset_name_list = []
        self._rf_name_list = []
        self._rf_nset_name_list = []
        self._rf_nsets = []
        self._load_nsets = [] # Nset definition of loads.
        self._load = self.generateLoadSetting() 

        # Assembly definition. 
        self._assembly_name = "assembly-1"
        self._instance_name = "instance-1"
        self._assembly_initial = ["*Assembly, name={}".format(self._assembly_name)] # Total list of Assembly definition. 
        self._instance = ["*Instance, name={}, part={}".format(self._instance_name, self._part_name),
                          "*End Instance"]

        self._ref_nodes_list = []
        
        self._fix_nset_name = "fix"
        self._fix_indices_list = fix_indices_list
        self._fix_nset = self.generateNset(self._fix_indices_list, self._fix_nset_name, self._instance_name) # Nset definition of fix BC. 
        self._loads_posi_indices_list = self._generateLoadPositions(self.loads_num, self._fix_indices_list) # Generate load positions. Randomly. For fixed mode: style="fix", input_posi_indices_list=[415, 470, 107]. 
        self._laplacian_initial_loads_posi = None # List. Containing the original position of concentrated forces.  
        self._laplacian_force_field = None # 2D Array of floats. Size: nSurfI * 3. The force field on the outer surface. 
        self._user_prescribed_force_field = user_prescribed_force_field # List of floats. Size: nSurfI * 3. The prescribed force field on the outer surface. Default: []. 

        self._surface_list = []
        self._coupling_list = [] 

        self._nset_boundary = [] # All nsets definitions in assembly. Boundary conditions
        self._assembly_end = ["*End Assembly"]
        self._assembly = self.generateAssembly()

        # Material. 
        self.material_type = material_type # String. Indicate material type. "linear"/"neo_hookean_fitting"/"neo_hookean_solid". 
        self._material_def_file_name = "" # Default: "". If there is a file of stress strain definition, please specify here (must not be ""). 
        self._material = self.generateMaterial(self.material_type)
        
        # Boundary condition. 
        self._boundary_initial = ["*Boundary"]
        self._boundary = self.generateBoundaryCondition_fixAll()

        # Step settings.  
        self.freq = int(self.autoIncrementNum / self.frameNum) # Int. The data frame extraction frequency (also refers to the number of increments. Extract one frame per "self.freq" increments). Especially for StaticSolver case.  
        
        self._step = ["*Step, name=step-1, nlgeom=YES, inc={}".format(self.autoIncrementNum),
                      "*Static",
                      "{}, {}, {}, {}".format(self.initIncrem, self.totalTime, 
                                              self.minIncrem, self.maxIncrem)] # Auto solver. 
        self._step_end = ["*End Step"]

        # Rest settings. 
        self._restart = ["*Restart, write, frequency=0"]
        self._output = ["*Output, field, variable=PRESELECT",
                        "*Output, history, variable=PRESELECT"]
        self._fil = ["*FILE FORMAT, ASCII",
                     "*node file, frequency={}".format(self.freq),
                     "U, COORD",
                     "*El file, frequency={}".format(self.freq),
                     "S, COORD"]
        self._resSettings = self._restart + self._output + self._fil
                
    
    def readFile(self, read_path):
        """
        Read files from specific path. 

        Parameters:
        ----------
            read_path: String. 
                Path of the original inp file.

        Return:
        ----------
            lines: List of strings. 
                The list of lines from the file. 
        """
        
        with open(read_path, "rt") as f: lines = f.read().splitlines()
        
        return lines
    

    def writeFile(self, write_status):
        """
        Write 'self.write_lines' into a new inp file. 

        Parameters:
        ----------
            write_status: String. 
                "Normal" / "Fast". 
                    "Normal": generate all definitions; 
                    "Fast": generate nodes and elements definition only. 
        """
        
        if write_status == "Normal":
            self._inputFile_lines_total = (self._header + self._part + self._assembly + 
                                        self._material + self._boundary + self._step + 
                                        self._load + self._resSettings + self._step_end)

            content = '\n'.join(self._inputFile_lines_total)
            
            with open(self.writePath, 'w') as f: f.write(content)
        
        elif write_status == "Fast":
            self._inputFile_lines_total = self._header + self._part

            content = '\n'.join(self._inputFile_lines_total)
            
            with open(self.writePath, 'w') as f: f.write(content)
        
        else:
            self.writeFile("Normal")
    

    def generatePart(self):
        """
        Generate part definition.  

        Returns:
        ----------
        The list collection of all sub-definition lists, including:
            part_initial: header part of "Part definition". 
            node: Node definition.
            elem: Element definition.
            elset_all: The elset containing all elements. For material definition specifically. 
            section: Section definition. 
            part_end: The endline of "Part definition". 
        """

        self.generateNodes(self.data_mat[self._node_variable_name], self._node)
        self.generateElements(self.data_mat[self._elem_variable_name], self._elem)

        self.nonlinearization()

        # Generate all element elset. 
        allElem_list, allElem_list_name = [], "allElems"
        for i in range(len(self._elem[1:])): allElem_list.append(str(i+1))
        self._elset_all = self.generateElset(allElem_list, allElem_list_name)

        # Generate Section. 
        self._section = self.generateSection(allElem_list_name, self._material_name)

        # Collection. 
        return (self._part_initial + self._node + self._elem + self._elset_all + 
                      self._section + self._part_end)


    def generateNodes(self, node_mat, target_node_list, specified_indices_list=[]):
        """
        Generate nodes information. 

        Parameters:
        ----------
            node_mat: 2D Array of ints. 
                The matrix containing the coordinates of the nodes to-be-defined under "*Node". 
            targer_node_list: List of strings. 
                The definition of node list. 
            specified_indices_list (optional): List of ints.
                List the indices of the input node list, following the exact order of the node_mat.  
                Default: []. 
        """

        for i in range(node_mat.shape[0]):
            if specified_indices_list == []: node_list_temp = ["{}".format(i+1)]
            else: node_list_temp = ["{}".format(specified_indices_list[i])]
            node_list_temp += [str(coord) for coord in list(node_mat[i,:])]
            target_node_list.append(', '.join(node_list_temp))
    
    
    def _extractOuterSurfaceNodes(self, faces_def_matrix, outer_surface_regionNum):
        """
        Extract the nodes on the outer surface of the geometry (for force application in next step). 

        Parameters:
        ----------
            faces_def_matrix: 2D Array of ints. 
                The definition of all faces, including the information of surface region number. 
            outer_surface_regionNum: Int. 
                The region number of outer surface of the geometry. 
        
        Returns:
        ----------
            outer_surface_nodes_list: List of ints. 
                The indices of nodes on the outer surface. Indexed from 1. Sorted. 
        """

        outer_surface_nodes_list = []

        for i in range(faces_def_matrix.shape[0]):
            if faces_def_matrix[i,0] == outer_surface_regionNum: # The region number of outer surface. 
                outer_surface_nodes_list += [int(ind) for ind in faces_def_matrix[i,1:]] # Indexed from 1. 
        
        outer_surface_nodes_list = list(set(outer_surface_nodes_list))
        outer_surface_nodes_list.sort()

        return outer_surface_nodes_list
    
    
    def generateElements(self, elem_mat, target_elem_list, specified_indices_list=[]):
        """
        Generate elements information. 

        Parameters:
        ----------
            elem_mat: 2D Array of ints. 
                The matrix containing the indices of each element to-be-defined under "*Element". 
            targer_elem_list: List of strings. 
                The definition of element list. 
            specified_indices_list (optional): List of ints.
                List the indices of the input element list, following the exact order of the elem_mat.  
                Default: []. 
        """

        for i in range(elem_mat.shape[0]):
            if specified_indices_list == []: elem_list_temp = ["{}".format(i+1)]
            else: elem_list_temp = ["{}".format(specified_indices_list[i])]

            elem_line_temp = [str(ind) for ind in list(elem_mat[i,:])]
            
            # Make sure the order of nodes for tetrahedron definition is counter-clockwise, otherwise resulting in negative volume. 
            ind_temp = elem_line_temp[1]
            elem_line_temp[1] = elem_line_temp[2]
            elem_line_temp[2] = ind_temp

            elem_list_temp += elem_line_temp
            target_elem_list.append(', '.join(elem_list_temp))
    

    def generateNset(self, node_list, nset_name, instance_name=None):
        """
        Generate node set information. 

        Parameters:
        ----------
        node_list: List of ints. 
            The list of nodes to be contained in the node list. 
        nset_name: String. 
            The name of the to-be-defined node list. 
        instance_name (optional): String. 
            The name of specified instance. 
            Only use in assembly definition. 
            Default: None. (Part cases)

        Returns:
        ----------
        nset: List of strings. 
            The definition of a specific nset. 

        """

        if instance_name == None: nset = ["*Nset, nset={}".format(nset_name)]
        else: nset = ["*Nset, nset={}, instance={}".format(nset_name, instance_name)]

        nset_line_temp, nset_string_temp = [], None

        for i, ind in enumerate(node_list):
            nset_line_temp.append(str(ind))

            if (i+1) % 10 == 0:
                nset_string_temp = ', '.join(nset_line_temp)
                nset.append(copy.deepcopy(nset_string_temp))
                nset_line_temp, nset_string_temp = [], None
        
        nset_string_temp = ', '.join(nset_line_temp)
        nset.append(copy.deepcopy(nset_string_temp))

        return nset
    
    
    def generateElset(self, elem_list, elset_name, instance_name=None):
        """
        Generate element set information.

        Parameters:
        ----------
        elem_list: List of ints. 
            The list of elements to be contained in the element list. 
        elset_name: String. 
            The name of the to-be-defined element list. 
        instance_name (optional): String. 
            The name of specified instance. 
            Only use in assembly definition. 
            Default: None. (Part cases)

        Returns:
        ----------
        elset: List of strings. 
            The definition of a specific elset.         
        """

        if instance_name == None: elset = ["*Elset, elset={}".format(elset_name)]
        else: elset = ["*Elset, elset={}, instance={}".format(elset_name, instance_name)]

        elset_line_temp, elset_string_temp = [], None

        for i, ind in enumerate(elem_list):
            elset_line_temp.append(str(ind))

            if (i+1) % 10 == 0:
                elset_string_temp = ', '.join(elset_line_temp)
                elset.append(copy.deepcopy(elset_string_temp))
                elset_line_temp, elset_string_temp = [], None

        elset_string_temp = ', '.join(elset_line_temp)
        elset.append(copy.deepcopy(elset_string_temp))

        return elset
    
    
    def generateSection(self, elset_name, material_name):
        """
        Generate section information. 

        Parameters:
        ----------
        elset_name: String. 
            The name of the elset to be assigned a section. 
        material_name: String. 
            The name of defined material. 

        Returns:
        ----------
        section: List of strings. 
            The definition of section. 
        """

        section = ["*Solid Section, elset={}, material={}".format(elset_name, material_name),
                   ","]

        return section
    
    
    def generateMaterial(self, material_type):
        """
        Generate lines for material definition. 

        Parameters:
        ----------
            material_type: String. 
                Indicate what type of material is used. 

        Returns:
        ----------
            material_lines: List of lines. 
                The lines of material definition. 
        """

        material_lines = ["*Material, name={}".format(self._material_name)]

        if material_type == "neo_hookean_fitting":
            stress_strain_lines = self._generateNeoHookeanFitting(self._modulus, (-0.3, 0.3), file_name=self._material_def_file_name)
            material_lines += ["*Hyperelastic, neo hooke, test data input, poisson={}".format(self._poisson_ratio), 
                               "*Uniaxial Test Data"]
            material_lines += stress_strain_lines
        
        elif material_type == "neo_hookean_solid":
            c10 = self._modulus / (4 * (1 + self._poisson_ratio))
            d1 = 6 * (1 - 2 * self._poisson_ratio) / self._modulus
            material_lines += ["*Hyperelastic, neo hooke",
                               "{}, {}".format(c10, d1)]

        elif material_type == "linear":
            material_lines += ["*Elastic",
                               "{}, {}".format(self._modulus, self._poisson_ratio)]
        
        else: material_lines = self.generateMaterial("linear")

        return material_lines


    def _generateNeoHookeanFitting(self, modulus, strain_range, file_name=""):
        """
        Import/Generate stress strain data for neo-Hookean material fitting. 

        Parameters:
        ----------
            modulus: Float. 
                The elastic modulus of material. 
            strain_range: Tuple of floats.
                Range for strain interpolation. 
            file_name (optional): String. 
                The name of stress strain data definition file.
                Default: "".  
        
        Returns:
        ----------
            stress_strain_lines: List of strings. 
                The lines of stress strain data. 
        """

        if file_name != "": return self.readFile(file_name)
        else:
            """
            Assumptions of neo-Hookean formulation:
                Incompressible (Poisson's ratio = ~0.5, small deformation).
                Undergoing uniaxial loading. 
                Formulation: sigma = 2*C*(stretch - 1/(stretch^2)). 
                E = 6*C. 
            """

            strain_data = np.linspace(strain_range[0], strain_range[1], 100)

            stretch_data = strain_data + 1.0
            stress_data = (self._modulus / 3.0) * (stretch_data - 1.0 / stretch_data**2) # Formulation. 

            stress_strain_lines = []

            for i in range(len(stress_data)):
                stress_strain_lines.append("%.6f, %.6f" % (stress_data[i], strain_data[i]))
            
            return stress_strain_lines

    
    def _generateLoadPositions(self, loads_num, fix_indices_list, style="random", input_posi_indices_list=[]):
        """
        Randomly generate positions of the load. 

        Parameters: 
        ----------
            loads_num: Int. 
                Number of loads.
            fix_indices_list: List of ints. 
                Indices of fixed nodes. 
            style (optional): String. 
                Indicate how to generate initial load positions. 
                "random" / "fix":
                    "random": Randomly generate load positions. 
                    "fix": Use the user input of initial load position indices.
                Default: "random". 
            input_posi_indices_list (optional): List of ints. 
                User input of initial load positions indices list. 
                Indexed from 1. 
                Default: []. 

        Returns:
        ----------
            loads_posi_indices_list: List of ints.
                Picked indices for load application positions. 
        """

        if style == "random":
            loads_posi_indices_list = []

            for i in range(loads_num):
                while(True):
                    load_posi_index_temp = random.choice(self._outer_surface_nodes_list) # Randomly chosen an outer surface node to apply load F(x, y, z). Indexed from 1. 
                    if load_posi_index_temp not in fix_indices_list: break # The randomly generated index cannot be one of the fixed nodes. 
                
                loads_posi_indices_list.append(load_posi_index_temp)

            return loads_posi_indices_list
        
        elif style == "fix": return input_posi_indices_list

        else: return self._generateLoadPositions(loads_num, fix_indices_list)
    

    def _generateLoadValues(self, output_dimension, load_scale, sampling_style="uniform"):
        """
        Randomly generate force values for load component definition.
        Using function: numpy.random.rand(). 

        Parameters:
        ----------
            output_dimension: Tuple of ints.
                The shape of output random array.  
                Size: 2*1. (dim1, dim2). 
            load_scale: Tuple of floats.
                Size: 2*1. (min_laod, max_laod) / (mean, deviation). 
            sampling_style (optional): String. 
                Indicating the type of sampling. 
                    "uniform": uniform distribution. 
                    "gaussian": Gaussian distribution. 
                Default: "uniform". 
        
        Returns:
        ----------
            load_result: Array of floats.
                Size: output_dimension. 
        """
        
        if sampling_style == "uniform":
            load_result = (np.random.rand(output_dimension[0], output_dimension[1]) * 2 - 1) * abs(load_scale[1] - load_scale[0])
            load_result = load_result.reshape(-1,1)

            for index, load_value_temp in enumerate(load_result):
                if load_value_temp < 0: load_result[index] -= self._load_scale[0]
                else: load_result[index] += self._load_scale[0]
            
            load_result = load_result.reshape(output_dimension[0], output_dimension[1])
        
        elif sampling_style == "gaussian":
            mean, deviation = load_scale[0], load_scale[1]
            load_result = np.random.normal(mean, deviation, size=output_dimension)
            load_result = load_result.reshape(-1,1)

            for index, load_value_temp in enumerate(load_result):
                if np.random.rand() <= 0.5: load_result[index] *= -1
            
            load_result = load_result.reshape(output_dimension[0], output_dimension[1])
        
        else: load_result = self._generateLoadValues(output_dimension, load_scale)

        return load_result

    
    def generateAssembly(self):
        """
        Generate assembly definition. 

        Returns:
        ----------
        The list collection of all sub-definition lists, including:
            assenbly_initial: Header of the assembly definition. 
            instance: The instance definition. 
            nset_boundary: The definition of BC related node set. 
            asssenbly_end: The endline of assembly definition. 
        """

        # Generate "self.loads_num" nsets, each of which has 1 node. 
                
        if self._isCoupleOn:
            for i, load_posi_index_temp in enumerate(self._loads_posi_indices_list): 
                ref_name_temp = "rf-{}".format(i+1)
                ref_nset_name_temp = "rf-{}-nset".format(i+1)
                self._rf_name_list.append(ref_name_temp)
                self._rf_nset_name_list.append(ref_nset_name_temp)

                # Generate assembly node definitions for reference points. 
                ref_node_list_temp = ["*Node"]
                ref_pt_coord_list_temp = [float(item) for item in self._node[load_posi_index_temp].split(',')[1:]]
                self.generateNodes(np.array(ref_pt_coord_list_temp).astype(float).reshape(1,-1), ref_node_list_temp, 
                                specified_indices_list=[i+1])
                self._ref_nodes_list += copy.deepcopy(ref_node_list_temp)
                rf_nset_list_temp = self._findCouplingNodes(load_posi_index_temp, self._coupling_neighbor_layers)

                # Generate reference point node sets. 
                self._load_nsets += self.generateNset([i+1], ref_name_temp)
                
                # Generate coupling constraint node sets. 
                self._rf_nsets += self.generateNset(rf_nset_list_temp, ref_nset_name_temp, 
                                                    self._instance_name)
            
            self.generateCoupling()
        
        else:
            if self._isLaplacianSmoothingOn:
                force_vector_temp = np.zeros(shape=(3*self._surface_nodes_num, 1))

                self._laplacian_initial_loads_posi = copy.deepcopy(self._loads_posi_indices_list)
                
                if self._initial_force_component_vector == []:
                    for load_posi_index_temp in self._loads_posi_indices_list:
                        force_vector_temp[(load_posi_index_temp-1)*3:load_posi_index_temp*3,:] = self._generateLoadValues((3,1), self._load_params_tuple, 
                                                                                                                        sampling_style=self._load_sampling_style)
                else: 
                    for load_posi_index_temp in self._loads_posi_indices_list:
                        force_vector_temp[(load_posi_index_temp-1)*3:load_posi_index_temp*3,:] = np.array(self._initial_force_component_vector).astype(float).reshape(3,1)
                
                laplacian_matrix, mass_matrix = self.data_mat[self._laplacian_variable_name], self.data_mat[self._massMatrix_variable_name]
                laplacian_matrix = self._laplacianMatrixShrink(laplacian_matrix, self._surface_nodes, self.data_mat["faces"], self._outer_surface_regionNum)

                force_vector_new = self._laplacianSmoothing(force_vector_temp, laplacian_matrix, mass_matrix, iter_num=self._laplacian_iter_num, 
                                                            smoothing_rate=self._smoothing_rate, laplacian_force_field=self._user_prescribed_force_field) # Size: (nSurfI x 3)*1. Fix force value: initial_BC_state="fix" (not recommended). 

                self._laplacian_force_field = force_vector_new.reshape(-1,3)

                self._loads_posi_indices_list = copy.deepcopy([(list(force_vector_new).index(item)//3)+1 for item in list(force_vector_new) if item != 0]) # Indexed from 1.
                self._loads_posi_indices_list = list(set(self._loads_posi_indices_list))
                self._loads_posi_indices_list.sort()

                for i, load_posi_index_temp in enumerate(self._loads_posi_indices_list):
                    load_nset_name_temp = "Load-{}".format(i+1)
                    self._loads_nset_name_list.append(load_nset_name_temp)
                    
                    self._load_nsets += self.generateNset([load_posi_index_temp], load_nset_name_temp, self._instance_name)
                
                self._load_nsets += self.generateNset(self._laplacian_initial_loads_posi, "Orig_loads_posi", self._instance_name)
                
                self._load = self.generateLoadSetting(force_list=list(force_vector_new.reshape(-1,1)))
            
            else:
                for i, load_posi_index_temp in enumerate(self._loads_posi_indices_list):
                    load_nset_name_temp = "Load-{}".format(i+1)
                    self._loads_nset_name_list.append(load_nset_name_temp)
                    
                    self._load_nsets += self.generateNset([load_posi_index_temp], load_nset_name_temp, self._instance_name)
     
        # Concatenate assembly subparts. 
        self._nset_boundary = self._nset_boundary + self._load_nsets + self._rf_nsets + self._fix_nset + self._surface_list + self._coupling_list

        return (self._assembly_initial + self._instance + self._ref_nodes_list + self._nset_boundary + self._assembly_end)


    def generateCoupling(self):
        """
        Generate coupling constriants for concentrated forces application. 
        """

        for index, rf_name in enumerate(self._rf_nset_name_list):
            self._surface_list += ["*Surface, type=NODE, name={}_CNS_, internal".format(rf_name),
                                   "{}, 1.".format(rf_name)]
            self._coupling_list += ["*Coupling, constraint name={}, ref node={}, surface={}_CNS_".format(self._rf_name_list[index], 
                                                                                                         self._rf_name_list[index],
                                                                                                         rf_name),
                                    "*{}".format(self._coupling_type)]
    

    def _findCouplingNodes(self, rf_node_ind, neighbor_layers):
        """
        Find the immediate neighbors of each specified node index. 

        Parameters:
        ----------
            rf_node_ind: Int. 
                The index of target node. 

        Returns:
        ----------
            rf_nset_list: List of ints (duplicated items removed). 
                "rf_node_ind"'s corresponding immediate neighbor nodes set. 
        """

        rf_nset_list, new_nodes_list, searched_nodes_list = [rf_node_ind], [rf_node_ind], []

        for j in range(neighbor_layers):
            for ind_temp in new_nodes_list:
                for i in range(len(self._triangle_nodes_list)):
                    if ind_temp in self._triangle_nodes_list[i]: 
                        rf_nset_list += copy.deepcopy(self._triangle_nodes_list[i])
                    else: continue
            
            searched_nodes_list += copy.deepcopy(new_nodes_list)
            rf_nset_list = list(set(copy.deepcopy(rf_nset_list)))
            new_nodes_list = [ind for ind in rf_nset_list if ind not in searched_nodes_list]

        # Avoid assigning same nodes to different coupled node sets. 
        for ind in rf_nset_list:
            if ind in self._coupled_list: rf_nset_list.remove(ind)
            else: self._coupled_list.append(ind)

        return rf_nset_list
    

    def generateBoundaryCondition_fixAll(self):
        """
        Generate fix boundary condition. 

        Returns:
        ----------
        The list collection of all sub-definition lists, including:
            boundary_initial: Header of boundary condition definition. 
            BC_list_temp: The detailed BC definition of boundary conditions. 
        """

        BC_list_temp = []
        for i in range(6): # 6: 6 DOFs (disp. + rot.); 3: 3 DOFs (disp.). 
            BC_list_temp.append("{}, {}, {}".format(self._fix_nset_name, i+1, i+1))
        
        return (self._boundary_initial + BC_list_temp)
    

    def generateLoadSetting(self, force_list=[]):
        """
        Generate load information. 

        Returns:
        ----------
            load_list: List of strings. 
                Definition of concentrated forces. 
            force_list (optional): List of forces (floats). 
                Size: loads_num * 3. 
                Default: [].
        """

        load_list = []

        if force_list == []: 
            force_list = list(self._generateLoadValues((self.loads_num*3, 1), self._load_params_tuple, sampling_style=self._load_sampling_style))
            
        force_list = np.array(force_list).astype(float).reshape(-1,3) # 2D Array of floats. Size: self._loads_num * 3. 

        if self._isCoupleOn:
            for j, rf_name in enumerate(self._rf_name_list): # Length: self._loads_num
                load_temp = ["*Cload, op=NEW"]

                for i in range(force_list.shape[1]): # 3: Three directions. 
                    load_temp.append("{}, {}, {}".format(rf_name, i+1, force_list[j,i]))
                
                load_list += copy.deepcopy(load_temp)

        else:
            for j, load_name in enumerate(self._loads_nset_name_list): # Length: length of self._loads_nset_name_list. 
                load_temp = ["*Cload"]

                for i in range(force_list.shape[1]): # 3: Three directions.
                    load_temp.append("{}, {}, {}".format(load_name, i+1, force_list[self._loads_posi_indices_list[j]-1,i]))
            
                load_list += copy.deepcopy(load_temp)
       
        return load_list
    
    
    def _laplacianMatrixShrink(self, laplacian_matrix, surface_nodes_list, faces_def_matrix, outer_surface_regionNum):
        """
        Assign zeros to the DOFs without force value applied. 

        Parameters:
        ----------
            laplacian_matrix: 2D Array of floats. 
                The surface's Laplacian for force smoothing. 
                Size: nSurfI*3 x nSurfI*3. 
            surface_nodes_list: List of ints. 
                All indices of nodes on all surfaces.
            faces_def_matrix: 2D Array of ints. 
                The definition of all faces, including the information of surface region number. 
            outer_surface_regionNum: Int. 
                The region number of outer surface of the geometry. 
        
        Returns:
        ----------
            laplacian_matrix: 2D Array of floats. 
                Laplacian with zeros assigned to the nodes not on the outer surfaces. 
                Size: nSurfI*3 x nSurfI*3. 
        """

        surface_nodes_list = [ind for ind in surface_nodes_list]
        
        outer_surface_nodes_list = self._extractOuterSurfaceNodes(faces_def_matrix, outer_surface_regionNum)
        other_surface_nodes_list = [ind for ind in surface_nodes_list if ind not in outer_surface_nodes_list]
        other_surface_nodes_list.sort()

        for ind in other_surface_nodes_list:
            laplacian_matrix[surface_nodes_list.index(ind)*3:(surface_nodes_list.index(ind)+1)*3,:] = 0.0
            laplacian_matrix[:,surface_nodes_list.index(ind)*3:(surface_nodes_list.index(ind)+1)*3] = 0.0

        return laplacian_matrix

    
    def _laplacianSmoothing(self, force_vector, laplacian_matrix, mass_matrix, iter_num=3, smoothing_rate=1e-4, initial_BC_state="", laplacian_force_field=[]):
        """
        Implement laplacian smoothing based on pre-calculated Laplacian matrix.
        Formulation: Forward Euler.
            F_(n+1) = (I + lambda*massMatrix*Laplacian) * F_n 

        Parameters:
        ----------
            force_vector: 1D Array of floats.
                With concentrated force values applied at the specidied nodes. 
                Size: (self._surface_nodes_num x 3) * 1. 
            laplacian_matrix: 2D Array of floats.
                Size: (self._surface_nodes_num x 3) * (self._surface_nodes_num x 3).
            mass_matrix: 2D Array of floats.
                Diagonal matrix. 
                Size: (self._surface_nodes_num x 3) * (self._surface_nodes_num x 3).
            iter_num (optional): Int. 
                The number of smoothing iterations.
                Default: 3. 
            smoothing_rate (optional): float.
                The coefficient that control the step size of smoothing. 
                Default: 1e-4. 
            initial_BC_state (optional): String. 
                Indicating whether to "fix" or "decay" the original concentrated force value.
                Default: "". Indicating smoothing including the original forces. 
            laplacian_force_field (optional): List of floats.
                The user-prescribed vector of laplacian force field. 
                Size: self._surface_nodes_num x 3. 
                Default: []. 

        Returns:
        ----------
            force_vector_new: 1D Array of floats. 
                The laplacian-smoothed force vector.
                Size: (self._surface_nodes_num x 3) * 1. 
        """

        if laplacian_force_field == []:
            force_vector_new = copy.deepcopy(force_vector)
            for i in range(iter_num): 
                force_vector_new += smoothing_rate * (laplacian_matrix @ force_vector_new) # Without mass matrix. 
                # force_vector_new += smoothing_rate * (mass_matrix @ laplacian_matrix @ force_vector_new) # With mass matrix (NOT recommended). 

                if initial_BC_state == "fix":
                    for j, value in enumerate(force_vector):
                        if value != 0:
                            force_vector_new[j] = value

        else: force_vector_new = np.array(laplacian_force_field).astype(float).reshape(len(laplacian_force_field),1)

        return force_vector_new


    def _computeMidPoint(self, ind_1, ind_2):
        """
        Compute the mid-point of the edge. 

        Parameters:
        ----------
            ind_1: Int. 
                The first index of the node pair. Indexed from 1. 
            ind_2: Int. 
                The second index of the node pair. Indexed from 1.
        
        Returns:
        ----------
            ind_mid: Int. 
                The index of the self._node. Index from 1. 
        """

        key_string_temp_1, key_string_temp_2 = "{}_{}".format(ind_1, ind_2), "{}_{}".format(ind_2, ind_1)

        if key_string_temp_1 in self._new_node_dict.keys(): return self._new_node_dict[key_string_temp_1]
        elif key_string_temp_2 in self._new_node_dict.keys(): return self._new_node_dict[key_string_temp_2]
        
        else:
            coord_temp_1 = np.array(self._node[ind_1].split(',')[1:]).astype(float).reshape(1,-1)
            coord_temp_2 = np.array(self._node[ind_2].split(',')[1:]).astype(float).reshape(1,-1)

            coord_temp_mid = (coord_temp_1 + coord_temp_2) / 2.0
            coord_mid_list = [str(item) for item in list(coord_temp_mid[0])]

            self._node_num = len(self._node)
            new_node_def_list_temp = copy.deepcopy([str(self._node_num)])
            new_node_def_list_temp += copy.deepcopy(coord_mid_list)
            self._node.append(', '.join(new_node_def_list_temp))
            self._new_node_list.append(', '.join(new_node_def_list_temp))

            self._new_node_dict[key_string_temp_1] = self._node_num
            self._new_node_dict[key_string_temp_2] = self._node_num

            return self._node_num
    

    def insertNode(self):
        """
        Insert one node (at the mid-point) of each edge. 
        Create C3D10 element structure. 
        """

        for index, elem_def_string in enumerate(self._elem[1:]):
            elem_node_list_temp = [int(ind) for ind in elem_def_string.split(',')[1:]]
            
            # Obtain the mid-point index in order. Assume tetrahedral element (C3D4). 
            mid_pt_ind_5 = self._computeMidPoint(elem_node_list_temp[0], elem_node_list_temp[1])
            mid_pt_ind_6 = self._computeMidPoint(elem_node_list_temp[1], elem_node_list_temp[2])
            mid_pt_ind_7 = self._computeMidPoint(elem_node_list_temp[0], elem_node_list_temp[2])
            mid_pt_ind_8 = self._computeMidPoint(elem_node_list_temp[0], elem_node_list_temp[3])
            mid_pt_ind_9 = self._computeMidPoint(elem_node_list_temp[1], elem_node_list_temp[3])
            mid_pt_ind_10 = self._computeMidPoint(elem_node_list_temp[2], elem_node_list_temp[3])

            elem_new_def_list_temp = [str(mid_pt_ind_5), 
                                      str(mid_pt_ind_6),
                                      str(mid_pt_ind_7),
                                      str(mid_pt_ind_8),
                                      str(mid_pt_ind_9),
                                      str(mid_pt_ind_10)]
            
            # Redefine the new C3D10 element in order. 
            elem_def_list_temp = copy.deepcopy(elem_def_string.split(',')) + copy.deepcopy(elem_new_def_list_temp)
            elem_def_string_temp = ', '.join(elem_def_list_temp)

            self._elem[index+1] = copy.deepcopy(elem_def_string_temp)

    
    def _triangleNodesCollection(self):
        """
        Collect all the nodes on each triangle (surface). 
        Need to be implemented after "self.insertNode()". 
        """

        for i in range(self._surface_mat.shape[0]):
            tri_temp = self._surface_mat[i,:]

            # Assuming all triangles on the surface of geometry. 
            middle_pts_list_temp = [self._computeMidPoint(tri_temp[0], tri_temp[1]),
                                    self._computeMidPoint(tri_temp[0], tri_temp[2]),
                                    self._computeMidPoint(tri_temp[1], tri_temp[2])]

            triangle_nodes_list_temp = list(copy.deepcopy(tri_temp)) + copy.deepcopy(middle_pts_list_temp)
            self._triangle_nodes_list.append(copy.deepcopy(triangle_nodes_list_temp)) # List of lists of ints. 

    
    def nonlinearization(self):
        """
        Nonlinearize the linear tetrahedral (CST) element to quadratic tetrahedral element.
        """

        self._elem_num = len(self._elem) - 1
        self._orig_node_num = len(self._node) - 1

        self.insertNode()
        self._triangleNodesCollection()

        self._node_num = len(self._node) - 1


def saveLog(file_name_list, elapsed_time_list, write_status, data_file_name, 
            sample_num, fix_indices_list, loads_num, load_sampling_type, load_param_tuple, 
            material_type, modulus, poisson_ratio, isCoupleOn, isLaplacianSmoothingOn, 
            coupling_type="", coupling_neighbor_layer_num=1, 
            laplacian_iter_num=5, laplacian_smoothing_rate=1e-4, write_path="nonlinear_case_generation.log"):
    """
    Save the nonlinear cases generation results into .log file. 

    Parameters:
    ----------
        file_name_list: List of strings. 
            Names of generated files. 
        elapsed_time_list: List of floats. 
            Elapsed time of generation for each input file.
            In exact order. 
        write_status: String. 
            Indicating the type of input file generation. 
            "Normal" / "Fast": 
                "Normal": generate all definitions; 
                "Fast": generate nodes and elements definition only.
        data_file_name: String. 
            The name of modeling data file.
            Format: .mat
        sample_num: Int. 
            Number of generated input files. 
        fix_indices_list: List of ints.
            Indices of fixed points.  
            Indexed from 1. 
        loads_num: Int. 
            The number of concentrated forces.
        load_sampling_type: String.
            The distribution type for force sampling. 
            "uniform" / "gaussian": 
                "uniform": uniform distribution with specified (min, max) range.
                "gaussian": gaussian distribution with specified (mean, dev) parameters. 
        load_param_tuple: tuple of floats.
            Parameters of load sampling. 
            load_sampling_type specific. 
        material_type: String. 
            The type of material.
            "linear" / "neo_hookean_solid" / "neo_hookean_fitting":
                "linear": linear elastic material.
                "neo_hookean_solid": neo-Hookean solid following the stain energy formulation. 
                "neo_hookean_fitting": neo-Hookean solid following the strass-strain curved fitted from user-input strss-strain data. 
        modulus: Float. 
            Elastic modulus of the material.
        poisson_ratio: Float.
            Poisson's ratio of the material. 
        isCoupleOn: Boolean indicator. 
            True: using coupling constraint for local force distribution.
            False: not using coupling constraint.   
        isLaplacianSmoothingOn: Boolean indicator.
            True: using Laplacian-Beltrami operator matrix to smooth the force distribution.
            False: not using Laplacian smoothing. 
        coupling_type (optional): String.
            The type of coupling constraint. 
            Default: "".
        coupling_neighbor_layer_num (optional): Int. 
            The number of neighbor layers to which the local force distributing goes. 
            Default: 1. 
        laplacian_iter_num (optional): Int. 
            The number of iteration for laplacian smoothing. 
            Default: 5.
        laplacian_smoothing_rate (optional): Float. 
            The rate of Laplacian smoothing. 
            Default: 1e-4.
        write_path (optional): String. 
            The path of to-be-written file. 
            Default: "nonlinear_case_generation.log". 
    """

    if isCoupleOn: isCoupleOn_status = "On"
    else: isCoupleOn_status = "Off"

    if isLaplacianSmoothingOn: isLaplacianSmoothingOn_status = "On"
    else: isLaplacianSmoothingOn_status = "Off"

    content = ["Data_file_name: {}".format(data_file_name), 
               "Sample_num = {}".format(sample_num),
               "Fixed_indices_list (indexed from 1): {}".format(fix_indices_list),
               "Material type: {}".format(material_type),
               "Elastic modulus = {} Pa".format(modulus),
               "Poisson's ratio = {}".format(poisson_ratio), 
               "Loads_num = {}".format(loads_num)]

    if load_sampling_type == "uniform":
        content += ["Load sampling type: {}".format(load_sampling_type),
                    "Load sampling range (min, max): {} N".format(load_param_tuple)]
    elif load_sampling_type == "gaussian":
        content += ["Load sampling type: {}".format(load_sampling_type),
                    "Load sampling parameters (mean, dev): {} N".format(load_param_tuple)]
    else:
        load_sampling_type = "uniform"
        content += ["Load sampling type: {}".format(load_sampling_type),
                    "Load sampling range (min, max): {} N".format(load_param_tuple)]

    content += ["Coupling constraint status: {}".format(isCoupleOn_status),
                "Laplacian smoothing status: {}".format(isLaplacianSmoothingOn_status)]

    if isCoupleOn:
        content += ["Coupling type: {}".format(coupling_type),
                    "Coupling neighbor layer numbers: {}".format(coupling_neighbor_layer_num)]
    
    if isLaplacianSmoothingOn:
        content += ["Laplacian smoothing iteration numbers = {}".format(laplacian_iter_num), 
                    "Laplacian smoothing rate = {}".format(laplacian_smoothing_rate)]
    
    content += ["----------------------------------------------------------",
                "Input file\t\tExport status\tGeneration status\tElapsed time/s"]

    elapsed_time_total = 0
    for i, file_name in enumerate(file_name_list):
        data_string_temp = "{}\t\t{}\t\tCompleted\t".format(file_name, write_status) + "\t%.8f" % (elapsed_time_list[i])
        content.append(data_string_temp)
        elapsed_time_total += elapsed_time_list[i]
    
    content += ["----------------------------------------------------------",
                "Total elapsed time: {} s".format(elapsed_time_total)]
    content = '\n'.join(content)
    
    with open(write_path, 'w') as f: f.write(content)


def main():
    abaqus_default_directory = "C:/temp" # Default working directory of Abaqus. 
    inp_folder = "inp_files"
    sample_nums = 1500
    data_file_path = "data_aorta.mat"
    node_variable_name, elem_variable_name = "NodeI", "EleI"
    results_folder_path_stress, results_folder_path_coor = "stress", "coor"
    material_type = "neo_hookean_solid" # "linear" / "neo_hookean_fitting" / "neo_hookean_solid". 
    fix_indices_list = [1148, 1156, 1169] # Specify the node to fix. At least 3. Indexed from 1. 
    write_status = "Normal" # String. "Normal" / "Fast". "Normal": generate all definitions; "Fast": generate nodes and elements definition only. 

    # ================================== Force interpolation related variables ================================== #
    force_field_mat_name = "force_field_data.mat"
    force_interpolation_folder = "inp_interpolation"
    isPrescribedForceOn = True # Boolean indicator. True: use prescribed force field; False: no specified force field. Default: False. 
    force_type = "random" # String. The type of prescribed force field. "interpolated": interpolated force fields; "random": weighted-summed force fields. 
    eigen_num_force, force_scalar = 20, 0.4 # Float. The scalar of force fields controlling the force magnitude -> deformation magnitude of the tumor in nonlinear solver. Unit: N. 
    # =========================================================================================================== #

    if isPrescribedForceOn:
        """
        The pipeline of generating interpolated force fields:
            1. Run "nonlinearCasesCreation.py" with 'isPrescribedForceOn = False' firstly. 
            2. Run "forceInterpolation.py" in the same directory. 
            3. Set 'isPrescribedForceOn = True', set 'force_type = "interpolated", then run "nonlinearCasesCreation.py" again. 
                Get input files with "*_interpolated.inp" in the folder 'force_interpolation_folder'. 
            4. Set 'isPrescribedForceOn = True', set 'force_type = "random", then run "nonlinearCasesCreation.py" again. 
                Get input files with "*_random.inp" in the folder 'force_interpolation_folder'. 
        """

        force_fields = (scipy.io.loadmat(force_field_mat_name)["force_field_interpolated"] if force_type == "interpolated" else 
                        scipy.io.loadmat(force_field_mat_name)["force_field_random"]) # Size: nSurfI*3 x sampleNum. Concatenated as xyzxyz...
        sample_nums = force_fields.shape[1]


    # Generate input file for Abaqus. 
    file_name_list, elapsed_time_list, force_field_matrix = [], [], None

    for i in range(sample_nums):
        start_time = time.time()

        if isPrescribedForceOn:
            if not os.path.isdir(force_interpolation_folder): os.mkdir(force_interpolation_folder)

            file_name_temp = ("{}_interpolated.inp".format(str(i+20001)) if force_type == "interpolated" else
                              "{}_random.inp".format(str(i+20001)))
            write_path = os.path.join(force_interpolation_folder, file_name_temp)

            force_field_prescribed_list = list(force_fields[:,i])

            inputFile_temp = inputFileGenerator(data_file_path, write_path, material_type, 
                                                fix_indices_list, node_variable_name, elem_variable_name, 
                                                user_prescribed_force_field=force_field_prescribed_list)

        else: 
            if not os.path.isdir(inp_folder): os.mkdir(inp_folder)

            file_name_temp = "{}.inp".format(str(i+20001))
            write_path = os.path.join(inp_folder, file_name_temp)

            inputFile_temp = inputFileGenerator(data_file_path, write_path, material_type, 
                                                fix_indices_list, node_variable_name, elem_variable_name)

        inputFile_temp.writeFile(write_status)
        end_time = time.time()
        elapsed_time = end_time - start_time

        file_name_list.append(file_name_temp)
        elapsed_time_list.append(elapsed_time)

        if i == 0: force_field_matrix = inputFile_temp._laplacian_force_field.reshape(-1,1)
        else: force_field_matrix = np.hstack((force_field_matrix, inputFile_temp._laplacian_force_field.reshape(-1,1)))

        # ============================ For force visualization only (sample_nums = 1) ============================ # 
        # print(inputFile_temp._laplacian_initial_loads_posi)
        # force_field = {"force_field": inputFile_temp._laplacian_force_field}
        # scipy.io.savemat("force_field.mat", force_field)
        # ======================================================================================================== # 

        print("Input_file: ", file_name_temp, "| Status:", write_status, "| Generation: Completed | Time: %.4f s" % (elapsed_time))

    saveLog(file_name_list, elapsed_time_list, write_status, data_file_path, sample_nums, 
            fix_indices_list, inputFile_temp.loads_num, inputFile_temp._load_sampling_style, inputFile_temp._load_params_tuple, 
            material_type, inputFile_temp._modulus, inputFile_temp._poisson_ratio, 
            inputFile_temp._isCoupleOn, inputFile_temp._isLaplacianSmoothingOn, 
            coupling_type=inputFile_temp._coupling_type, coupling_neighbor_layer_num=inputFile_temp._coupling_neighbor_layers, 
            laplacian_iter_num=inputFile_temp._laplacian_iter_num, laplacian_smoothing_rate=inputFile_temp._smoothing_rate, 
            write_path="nonlinear_case_generation.log")  

    if not isPrescribedForceOn: weight_matrix = (2.0 * np.random.rand(eigen_num_force, 3*sample_nums) - 1.0) # Distinct random weights corresponding to each laplacian-force-field.
    else: weight_matrix = scipy.io.loadmat(force_field_mat_name)["weight_matrix"] # Distinct random force field for each laplacian-force-field.

    mdict = {"fix_indices_list": fix_indices_list,
             "orig_data_file_name": data_file_path,
             "orig_config_var_name": node_variable_name,
             "inp_folder": inp_folder if not isPrescribedForceOn else force_interpolation_folder, # The folder containing input files. 
             "current_directory": os.getcwd(),
             "results_folder_path_stress": results_folder_path_stress,
             "results_folder_path_coor": results_folder_path_coor,
             "original_node_number": inputFile_temp._orig_node_num,
             "couple_region_num": inputFile_temp._couple_region_num,
             "force_field_matrix": force_field_matrix, # The force field matrix of all generated samples. Size: nSurfI*3 x sampleNum_total. 
             "weight_matrix": weight_matrix, "force_scalar_coeff": force_scalar, # The randomly generated matrix for force fields' reconstruction. Size: eigen_num x (3*sample_num). 
             "eigen_number_force": eigen_num_force, # Int. The eigenmode number of force field reconstruction. (Used only in force field interpolation)
             "alpha_indexing_vector": np.zeros(shape=(sample_nums, 1)) if not isPrescribedForceOn else scipy.io.loadmat(force_field_mat_name)["alpha_indexing_vector"]
            }

    scipy.io.savemat("training_parameters_transfer.mat", mdict)
    
    # np.save(os.path.join(abaqus_default_directory, "training_parameters_transfer.npy"), mdict, fix_imports=True)

    # np.savez(os.path.join(abaqus_default_directory, "training_parameters_transfer.npz"), 
    #          fix_indices_list=fix_indices_list,
    #          orig_data_file_name=data_file_path,
    #          orig_config_var_name=node_variable_name,
    #          inp_folder=inp_folder,
    #          current_directory=os.getcwd(),
    #          results_folder_path_stress=results_folder_path_stress,
    #          results_folder_path_coor=results_folder_path_coor)


if __name__ == "__main__":
    main()

[(1187, 'time.time', 'time.time', ([], {}), False, 'import time\n'), (1212, 'time.time', 'time.time', ([], {}), False, 'import time\n'), (1244, 'os.getcwd', 'os.getcwd', ([], {}), False, 'import os\n'), (393, 'copy.deepcopy', 'copy.deepcopy', (['nset_string_temp'], {}), False, 'import copy\n'), (433, 'copy.deepcopy', 'copy.deepcopy', (['elset_string_temp'], {}), False, 'import copy\n'), (529, 'numpy.linspace', 'np.linspace', (['strain_range[0]', 'strain_range[1]', '(100)'], {}), True, 'import numpy as np\n'), (761, 'copy.deepcopy', 'copy.deepcopy', (['new_nodes_list'], {}), False, 'import copy\n'), (906, 'copy.deepcopy', 'copy.deepcopy', (['force_vector'], {}), False, 'import copy\n'), (990, 'copy.deepcopy', 'copy.deepcopy', (['elem_def_string_temp'], {}), False, 'import copy\n'), (1194, 'os.path.join', 'os.path.join', (['force_interpolation_folder', 'file_name_temp'], {}), False, 'import os\n'), (1206, 'os.path.join', 'os.path.join', (['inp_folder', 'file_name_temp'], {}), False, 'import os\n'), (1252, 'numpy.zeros', 'np.zeros', ([], {'shape': '(sample_nums, 1)'}), True, 'import numpy as np\n'), (622, 'numpy.random.normal', 'np.random.normal', (['mean', 'deviation'], {'size': 'output_dimension'}), True, 'import numpy as np\n'), (662, 'copy.deepcopy', 'copy.deepcopy', (['ref_node_list_temp'], {}), False, 'import copy\n'), (676, 'numpy.zeros', 'np.zeros', ([], {'shape': '(3 * self._surface_nodes_num, 1)'}), True, 'import numpy as np\n'), (678, 'copy.deepcopy', 'copy.deepcopy', (['self._loads_posi_indices_list'], {}), False, 'import copy\n'), (818, 'copy.deepcopy', 'copy.deepcopy', (['load_temp'], {}), False, 'import copy\n'), (827, 'copy.deepcopy', 'copy.deepcopy', (['load_temp'], {}), False, 'import copy\n'), (952, 'copy.deepcopy', 'copy.deepcopy', (['coord_mid_list'], {}), False, 'import copy\n'), (987, 'copy.deepcopy', 'copy.deepcopy', (['elem_new_def_list_temp'], {}), False, 'import copy\n'), (1007, 'copy.deepcopy', 'copy.deepcopy', (['middle_pts_list_temp'], {}), False, 'import copy\n'), (1008, 'copy.deepcopy', 'copy.deepcopy', (['triangle_nodes_list_temp'], {}), False, 'import copy\n'), (1190, 'os.path.isdir', 'os.path.isdir', (['force_interpolation_folder'], {}), False, 'import os\n'), (1190, 'os.mkdir', 'os.mkdir', (['force_interpolation_folder'], {}), False, 'import os\n'), (1203, 'os.path.isdir', 'os.path.isdir', (['inp_folder'], {}), False, 'import os\n'), (1203, 'os.mkdir', 'os.mkdir', (['inp_folder'], {}), False, 'import os\n'), (1237, 'numpy.random.rand', 'np.random.rand', (['eigen_num_force', '(3 * sample_nums)'], {}), True, 'import numpy as np\n'), (389, 'copy.deepcopy', 'copy.deepcopy', (['nset_string_temp'], {}), False, 'import copy\n'), (429, 'copy.deepcopy', 'copy.deepcopy', (['elset_string_temp'], {}), False, 'import copy\n'), (574, 'random.choice', 'random.choice', (['self._outer_surface_nodes_list'], {}), False, 'import random\n'), (762, 'copy.deepcopy', 'copy.deepcopy', (['rf_nset_list'], {}), False, 'import copy\n'), (1007, 'copy.deepcopy', 'copy.deepcopy', (['tri_temp'], {}), False, 'import copy\n'), (611, 'numpy.random.rand', 'np.random.rand', (['output_dimension[0]', 'output_dimension[1]'], {}), True, 'import numpy as np\n'), (626, 'numpy.random.rand', 'np.random.rand', ([], {}), True, 'import numpy as np\n'), (758, 'copy.deepcopy', 'copy.deepcopy', (['self._triangle_nodes_list[i]'], {}), False, 'import copy\n'), (809, 'numpy.array', 'np.array', (['force_list'], {}), True, 'import numpy as np\n'), (916, 'numpy.array', 'np.array', (['laplacian_force_field'], {}), True, 'import numpy as np\n'), (660, 'numpy.array', 'np.array', (['ref_pt_coord_list_temp'], {}), True, 'import numpy as np\n'), (686, 'numpy.array', 'np.array', (['self._initial_force_component_vector'], {}), True, 'import numpy as np\n')]