diff --git a/__init__.py b/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/__init__.pyc b/__init__.pyc new file mode 100644 index 0000000..02f92be Binary files /dev/null and b/__init__.pyc differ diff --git a/__pycache__/__init__.cpython-37.pyc b/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 0000000..a812b4d Binary files /dev/null and b/__pycache__/__init__.cpython-37.pyc differ diff --git a/__pycache__/analysis.cpython-37.pyc b/__pycache__/analysis.cpython-37.pyc new file mode 100644 index 0000000..bcf268d Binary files /dev/null and b/__pycache__/analysis.cpython-37.pyc differ diff --git a/__pycache__/cross_section.cpython-37.pyc b/__pycache__/cross_section.cpython-37.pyc new file mode 100644 index 0000000..3f9961d Binary files /dev/null and b/__pycache__/cross_section.cpython-37.pyc differ diff --git a/__pycache__/presg.cpython-37.pyc b/__pycache__/presg.cpython-37.pyc new file mode 100644 index 0000000..2bacc7d Binary files /dev/null and b/__pycache__/presg.cpython-37.pyc differ diff --git a/__pycache__/sg.cpython-37.pyc b/__pycache__/sg.cpython-37.pyc new file mode 100644 index 0000000..2987494 Binary files /dev/null and b/__pycache__/sg.cpython-37.pyc differ diff --git a/__pycache__/utils.cpython-37.pyc b/__pycache__/utils.cpython-37.pyc new file mode 100644 index 0000000..c692c68 Binary files /dev/null and b/__pycache__/utils.cpython-37.pyc differ diff --git a/analysis.py b/analysis.py new file mode 100644 index 0000000..620a15b --- /dev/null +++ b/analysis.py @@ -0,0 +1,124 @@ +import os +import sys +import traceback +import datetime as dt +import subprocess as sbp +import msgpi.presg as psg +import msgpi.sg as sgm +import msgpi.io.iosc as miosc +import msgpi.io.iovabs as miovabs + + +def solve(sg_xml, analysis, solver, scrnout=True): + """ + + Parameters + ---------- + sg_xml : str + File name of SG design parameters (XML format) + analysis : str + Analysis to be carried out + h - homogenization + d - dehomogenization/localization/recover + f - initial failure strength + fe - initial failure envelope + fi - initial failure indices and strength ratios + solver : str + Format of the generated input file ('vabs' or 'swiftcomp') + """ + + # Preprocess + sg_in, smdim = psg.preSG(sg_xml, analysis, solver) + + # Solve + run(sg_in, analysis, solver, smdim, scrnout) + + # Parse results + print(' - reading results...') + if analysis == 'h': + if solver == 'vabs': + sm = miovabs.readVABSOutHomo(sg_in + '.k') + elif solver == 'swiftcomp': + sm = miosc.readSCOutHomo(sg_in + '.k', smdim) + return sm + elif analysis == 'd': + pass + elif 'f' in analysis: + results = miosc.readSCOutFailure(sg_in, analysis) + return results + + +def run(input_name, analysis, solver, smdim, scrnout=True): + """ Run codes + + Parameters + ---------- + solver : str + Excect command string of the code ('vabs' or 'swiftcomp') + input_name : str + Name of the input file + analysis : str + Analysis to be carried out + h - homogenization + d or l - dehomogenization/localization/recover + f - initial failure strength + fe - initial failure envelope + fi - initial failure indices and strength ratios + smdim : int + Dimension of the macroscopic structural model + + :param scrnout: Print solver messages + :type scrnout: bool + """ + try: + analysis_long = { + 'h': 'homogenization', + 'ha': 'homogenization aperiodic', + 'd': 'recover', + 'l': 'dehomogenization', + 'la': 'dehomogenization aperiodic', + 'lg': 'dehomogenization gmsh format', + 'lag': 'dehomogenization aperiodic gmsh format', + 'f': 'initial failure strength', + 'fe': 'initial failure envelope', + 'fi': 'initial failure indices strength ratios' + } + msg = ' - running {cn} for {an}...\n' + + if solver == 'vabs': + solver = solver + 'iii' + cmd = [solver, input_name] + + if solver == 'swiftcomp': + cmd.append(str(smdim) + 'D') + if 'd' in analysis: + analysis = analysis.replace('d', 'l') + cmd.append(analysis.upper()) + + cmd = ' '.join(cmd) + + # try: + if scrnout: + sys.stdout.write(msg.format( + cn=solver, an=analysis_long[analysis] + )) + sys.stdout.flush() + sbp.call(cmd) + else: + FNULL = open(os.devnull, 'w') + sbp.call(cmd, stdout=FNULL, stderr=sbp.STDOUT) + # except: + # sys.stdout.write('failed\n') + # sys.stdout.flush() + # return + # else: + # sys.stdout.write('done\n') + # sys.stdout.flush() + except: + e = traceback.format_exc() + print(e) + + +def runBatch(): + + return diff --git a/analysis.pyc b/analysis.pyc new file mode 100644 index 0000000..18137af Binary files /dev/null and b/analysis.pyc differ diff --git a/cross_section.py b/cross_section.py new file mode 100644 index 0000000..ffdb080 --- /dev/null +++ b/cross_section.py @@ -0,0 +1,131 @@ +import csv +import os +import sys +import numpy as np +import xml.etree.ElementTree as et + + +class Layer(object): + def __init__(self): + self.material = '' + self.thickness = 0 + self.lamina = None + self.angle = 0 + self.begin = 0 + self.end = 0 + + def __str__(self): + return 'begin: {0}, end: {1}, material: {2}, thickness: {3}, angle: {4}'.format( + self.begin, self.end, self.material, self.thickness, self.angle + ) + + def readIXGENLine(self, line, cs, fmt=0, tb=''): + line = line.split() + if (fmt == 0): + self.material = ' '.join(line[5:-3]) + self.thickness = line[-3] + self.angle = float(line[-2]) + if (cs == 1): + self.begin = float(line[1]) + self.end = float(line[2]) + elif (cs == 2): + self.begin = float(line[3]) + self.end = float(line[4]) + elif (fmt == 1): + self.material = ' '.join(line[4:-1]) + self.thickness = line[-1] + if (cs == 1): + if (tb == 'top'): + self.end = float(line[0]) + elif (tb == 'bottom'): + self.end = float(line[1]) + elif (cs == 2): + if (tb == 'top'): + self.end = float(line[2]) + elif (tb == 'bottom'): + self.end = float(line[3]) + elif (fmt == 2): + self.material = ' '.join(line[1:-2]) + self.thickness = line[-2] + self.angle = float(line[-1]) + + +class Layup(object): + def __init__(self, name=''): + self.name = name + self.layers = [] + + +class CrossSection(object): + """ Stores all information of a cross section. + + This object is used to store and pass data of a cross section + between VABS/SC, GEBT, and preprocessors. + """ + + def __init__(self, name): + self.name = name + + self.name_template = '' #: Base name of the main input file template + self.name_baseline = '' #: Base name of the base lines file + self.name_layup = '' #: Base name of the layups file + self.fn_vabs_in = self.name + '.sg' #: File name of the VABS input + self.fn_sc_in = self.name + '.sg' #: File name of the SwiftComp input + # self.fn_gmsh_msh = self.name + '.msh' #: File name of the Gmsh mesh file + + self.num_layertypes = 0 #: Number of layer types + self.layertypes = {} #: Layer types {lid: [mid, ply angle], ...} + self.element_layertype = {} + self.element_theta1 = {} + + self.pitch = 0.0 #: Pitch angle of the cross section + self.k11 = 0.0 #: Initial twist k_{11} + self.k12 = 0.0 #: Initial curvature k_{12} + self.k13 = 0.0 #: Initial curvature k_{13} + self.cos11 = 1.0 #: Obliqueness. Cosine of the angle between y_1 and x_1 + self.cos21 = 0.0 #: Obliqueness. Cosine of the angle between y_2 and x_1 + + self.trapeze = 0 #: Flag for trapeze effect + self.vlasov = 0 #: Flag for Vlasov model + + # Effective properties + self.mc = None #: Mass center + self.gc = None #: Geometric center + self.tc = None #: Tension center + self.sc = None #: Shear center + + self.mass = None #: Mass matrix (6x6) + self.masc = None #: Mass matrix at mass center (6x6) + self.eden = None #: Effective density (mass per unit length) + self.ecsm = None #: Effective classical stiffness matrix (4x4) + self.eccm = None #: Effective classical compliance matrix (4x4) + self.etsm = None #: Effective timoshenko stiffness matrix (6x6) + self.etcm = None #: Effective timoshenko compliance matrix (6x6) + + # Global results for recovery/dehomogenization + self.gdisplacements = np.zeros(3) #: Global displacements [u1, u2, u3] + #: Global rotations [[c11, c12, c13], [c21, c22, c23], [c31, c32, c33]] + self.grotations = np.eye(3) + self.gforces = np.zeros(3) #: Global forces (VABS) [F1, F2, F3] + self.gmoments = np.zeros(3) #: Global moments (VABS) [M1, M2, M3] + #: Global distributed forces (VABS) [[f1, f2, f3], [d1f1, d1f2, d1f3], [d2f1, d2f2, d2f3], [d3f1, d3f2, d3f3]] + self.gdforces = np.zeros((4, 3)) + #: Global distributed moments (VABS) [[m1, m2, m3], [d1m1, d1m2, d1m3], [d2m1, d2m2, d2m3], [d3m1, d3m2, d3m3]] + self.gdmoments = np.zeros((4, 3)) + #: Indicate whether the generalized loads are stresses or strains (SwiftComp) + self.gmeasure = 'stress' + #: Global straints (SwiftComp) [] + self.gloads = np.zeros(6) + + #: Load for the horizontal axis for failure envelope (SwiftComp) + self.feaxis1 = '' + #: Load for the vertical axis for failure envelope (SwiftComp) + self.feaxis2 = '' + #: Number of evaluation points along each axis for failure envelope (SwiftComp) + self.fendiv = 10 + + + + + + diff --git a/io/__init__.py b/io/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/io/__init__.pyc b/io/__init__.pyc new file mode 100644 index 0000000..8e25192 Binary files /dev/null and b/io/__init__.pyc differ diff --git a/io/__pycache__/__init__.cpython-37.pyc b/io/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 0000000..7d2ca6b Binary files /dev/null and b/io/__pycache__/__init__.cpython-37.pyc differ diff --git a/io/__pycache__/__init__.cpython-38.pyc b/io/__pycache__/__init__.cpython-38.pyc new file mode 100644 index 0000000..bfdb5ba Binary files /dev/null and b/io/__pycache__/__init__.cpython-38.pyc differ diff --git a/io/__pycache__/iosc.cpython-37.pyc b/io/__pycache__/iosc.cpython-37.pyc new file mode 100644 index 0000000..9842070 Binary files /dev/null and b/io/__pycache__/iosc.cpython-37.pyc differ diff --git a/io/__pycache__/iovabs.cpython-37.pyc b/io/__pycache__/iovabs.cpython-37.pyc new file mode 100644 index 0000000..2abba63 Binary files /dev/null and b/io/__pycache__/iovabs.cpython-37.pyc differ diff --git a/io/__pycache__/utils.cpython-37.pyc b/io/__pycache__/utils.cpython-37.pyc new file mode 100644 index 0000000..ec27d9a Binary files /dev/null and b/io/__pycache__/utils.cpython-37.pyc differ diff --git a/io/iosc.py b/io/iosc.py new file mode 100644 index 0000000..9996dc7 --- /dev/null +++ b/io/iosc.py @@ -0,0 +1,639 @@ +import os +import numpy as np +import msgpi.io.utils as utl +import msgpi.sg as sgm + + +# ==================================================================== +# +# +# READERS +# +# +# ==================================================================== + +def readSCIn(fn_sg, smdim): + """ Read data from the SwiftComp input file + + :param fn_sg: File name of the SwiftComp input file + :type fn_sg: string + """ + + fn_base, fn_extn = os.path.splitext(fn_sg) + name = os.path.basename(fn_base) + sg = sgm.StructureGene(name, 3, smdim) + + count = 0 + count_node = 1 + count_element = 1 + count_layertype = 1 + count_material = 1 + line_material = 0 + extra_head = 0 + if smdim == 1: + extra_head = 3 + elif smdim == 2: + extra_head = 2 + head = 2 # at least 3 lines in the head (flag) part + with open(fn_sg, 'r') as fin: + for li, line in enumerate(fin): + line = line.strip() + if line == '\n' or line == '': + continue + + count += 1 + line = line.split() + if count <= (extra_head + head): + # Read head + if smdim == 1: + if count == 1: + # model + line = list(map(int, line)) + sg.model = line[0] + elif count == 2: + # initial twist/curvatures + line = list(map(float, line)) + sg.initial_twist = line[0] + sg.initial_curvature = [line[1], line[2]] + elif count == 3: + # oblique + line = list(map(float, line)) + sg.oblique = line + elif smdim == 2: + if count == 1: + # model + line = list(map(int, line)) + sg.model = line[0] + elif count == 2: + # initial twist/curvature + line = list(map(float, line)) + sg.initial_curvature = line + + if count == (extra_head + head - 1): + # analysis elem_flag trans_flag temp_flag + line = list(map(int, line)) + sg.analysis = line[0] + sg.degen_element = line[1] + sg.trans_element = line[2] + sg.nonuniform_temperature = line[3] + elif count == (extra_head + head): + # nsg nnode nelem nmate nslave nlayer + line = list(map(int, line)) + sg.sgdim = line[0] + num_nodes = line[1] + num_elements = line[2] + num_materials = line[3] + num_layertypes = line[4] + continue + elif count_node <= num_nodes: + # Read nodal coordinates + sg.nodes[int(line[0])] = list(map(float, line[1:])) + count_node += 1 + continue + elif count_element <= num_elements: + # Read element material/layer connectivities + line = list(map(int, line)) + eid = line[0] + lyrtp = line[1] + + sg.elementids.append(eid) + sg.elements[eid] = [i for i in line[2:] if i != 0] + sg.elementids2d.append(eid) + + if lyrtp not in sg.prop_elem.keys(): + sg.prop_elem[lyrtp] = [] + sg.prop_elem[lyrtp].append(eid) + sg.elem_prop[eid] = lyrtp + + count_element += 1 + continue + elif (sg.trans_element != 0) and (count_element <= (2 * num_elements)): + # Read element orientation (a b c) + eid = int(line[0]) + a = list(map(float, line[1:4])) + b = list(map(float, line[4:7])) + c = list(map(float, line[7:10])) + sg.elem_orient[eid] = [a, b, c] + + count_element += 1 + continue + elif count_layertype <= num_layertypes: + # Read layer types + ltid = int(line[0]) + mid = int(line[1]) + theta3 = float(line[2]) + sg.mocombos[ltid] = [mid, theta3] + count_layertype += 1 + continue + elif count_material <= num_materials: + # Read materials + if line_material == 0: + line = list(map(int, line)) + mid = line[0] + mtype = line[1] + num_temp = line[2] + sg.materials[mid] = sgm.MaterialSection() + sg.materials[mid].eff_props[3]['type'] = mtype + line_material += 1 + continue + elif line_material == 1: + line = list(map(float, line)) + sg.materials[mid].eff_props[3]['density'] = dens + line_material += 1 + continue + + if mtype == 0: + # isotropic + if line_material == 2: + e = float(line[0]) + nu = float(line[1]) + sg.materials[mid].eff_props[3]['constants']['e'] = e + sg.materials[mid].eff_props[3]['constants']['nu'] = nu + line_material = 99 + continue + elif mtype == 1: + # orthotropic + if line_material == 2: + # pe = [] + e = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['e1'] = e[0] + sg.materials[mid].eff_props[3]['constants']['e2'] = e[1] + sg.materials[mid].eff_props[3]['constants']['e3'] = e[2] + # for ei in e: + # pe.append(ei) + # sg.materials[mid].property_elastic += e + line_material += 1 + continue + elif line_material == 3: + g = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['g12'] = g[0] + sg.materials[mid].eff_props[3]['constants']['g13'] = g[1] + sg.materials[mid].eff_props[3]['constants']['g23'] = g[2] + # for gi in g: + # pe.append(gi) + # self.materials[mid].property_elastic += g + line_material += 1 + continue + elif line_material == 4: + nu = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['nu12'] = nu[0] + sg.materials[mid].eff_props[3]['constants']['nu13'] = nu[1] + sg.materials[mid].eff_props[3]['constants']['nu23'] = nu[2] + # for nui in nu: + # pe.append(nui) + # sg.materials[mid].property_elastic.append(pe) + line_material = 99 + continue + elif mtype == 2: + # anisotropic + continue + + if line_material == 99: + line_material = 0 + count_material += 1 # next material + continue + continue + else: + # omega + sg.omega = float(line[0]) + + return sg + + +def readSCOutHomo(fn, smdim): + """Read SwiftComp homogenization results. + + :param fn: SwiftComp output file (e.g. example.sg.k) + :type fn: string + + :param smdim: Dimension of the structural model + :type smdim: int + """ + linesRead = [] + keywordsIndex = {} + + with open(fn, 'r') as fin: + for lineNumber, line in enumerate(fin): + linesRead.append(line) + if 'The Effective Stiffness Matrix' in line: + keywordsIndex['The Effective Stiffness Matrix'] = lineNumber + if 'The Effective Compliance Matrix' in line: + keywordsIndex['The Effective Compliance Matrix'] = lineNumber + if smdim == 1: + if 'Timoshenko Stiffness' in line: + keywordsIndex['Timoshenko Stiffness'] = lineNumber + if 'Timoshenko Compliance' in line: + keywordsIndex['Timoshenko Compliance'] = lineNumber + if 'Shear Center Location' in line: + keywordsIndex['Shear Center Location'] = lineNumber + elif smdim == 2: + if 'In-Plane Properties' in line: + keywordsIndex['In-Plane Properties'] = lineNumber + if 'Flexural Properties' in line: + keywordsIndex['Flexural Properties'] = lineNumber + elif smdim == 3: + if 'The Engineering Constants' in line: + keywordsIndex['The Engineering Constants'] = lineNumber + if 'Effective Density' in line: + keywordsIndex['Effective Density'] = lineNumber + + sm = sgm.MaterialSection(smdim) + + # ems = 6 # effective matrix size + + if smdim == 1: + ems = 4 + if 'The Effective Stiffness Matrix' in keywordsIndex.keys(): + indexStiffness = keywordsIndex['The Effective Stiffness Matrix'] + sm.eff_props[1]['stiffness']['classical'] = utl.textToMatrix( + linesRead[indexStiffness + 2:indexStiffness + 2 + ems]) + # else: + # print('No effective stiffness matrix found.') + + if 'The Effective Compliance Matrix' in keywordsIndex.keys(): + indexCompliance = keywordsIndex['The Effective Compliance Matrix'] + sm.eff_props[1]['compliance']['classical'] = utl.textToMatrix( + linesRead[indexCompliance + 2:indexCompliance + 2 + ems]) + # else: + # print('No effective compliance matrix found.') + + if 'Timoshenko Stiffness' in keywordsIndex.keys(): + indexStiffness = keywordsIndex['Timoshenko Stiffness'] + sm.eff_props[1]['stiffness']['refined'] = utl.textToMatrix( + linesRead[indexStiffness + 2:indexStiffness + 8]) + # else: + # print('No effective timoshenko stiffness matrix found.') + + if 'Timoshenko Compliance' in keywordsIndex.keys(): + indexCompliance = keywordsIndex['Timoshenko Compliance'] + sm.eff_props[1]['compliance']['refined'] = utl.textToMatrix( + linesRead[indexCompliance + 2:indexCompliance + 8]) + # else: + # print('No effective timoshenko compliance matrix found.') + + if 'Shear Center Location' in keywordsIndex.keys(): + index_sc = keywordsIndex['Shear Center Location'] + sm.eff_props[1]['shearcenter'] = list(map( + float, linesRead[index_sc+2].strip().split() + )) + # else: + # print('No shear center location found.') + + line = linesRead[keywordsIndex['Effective Density']] + line = line.strip().split('=') + sm.eff_props[1]['density'] = float(line[-1].strip()) + + elif smdim == 2: + ems = 6 + + try: + indexStiffness = keywordsIndex['The Effective Stiffness Matrix'] + sm.eff_props[2]['stiffness']['classical'] = utl.textToMatrix( + linesRead[indexStiffness + 2:indexStiffness + 2 + ems]) + except KeyError: + print('No effective stiffness matrix found.') + + try: + indexCompliance = keywordsIndex['The Effective Compliance Matrix'] + sm.eff_props[2]['compliance']['classical'] = utl.textToMatrix( + linesRead[indexCompliance + 2:indexCompliance + 2 + ems]) + except KeyError: + print('No effective compliance matrix found.') + + try: + index = keywordsIndex['In-Plane Properties'] + for line in linesRead[index + 2:index + 8]: + line = line.strip() + line = line.split('=') + label = line[0].strip().lower() + value = float(line[-1].strip()) + sm.eff_props[2]['constants']['inplane'][label] = value + except KeyError: + print('No in-plane properties found.') + + try: + index = keywordsIndex['Flexural Properties'] + for line in linesRead[index + 2:index + 8]: + line = line.strip() + line = line.split('=') + label = line[0].strip().lower() + value = float(line[-1].strip()) + sm.eff_props[2]['constants']['flexural'][label] = value + except KeyError: + print('No flexural properties found.') + + line = linesRead[keywordsIndex['Effective Density']] + line = line.strip().split('=') + sm.eff_props[2]['density'] = float(line[-1].strip()) + + elif smdim == 3: + ems = 6 + + try: + indexStiffness = keywordsIndex['The Effective Stiffness Matrix'] + sm.eff_props[3]['stiffness'] = utl.textToMatrix( + linesRead[indexStiffness + 2:indexStiffness + 2 + ems]) + except KeyError: + print('No effective stiffness matrix found.') + + try: + indexCompliance = keywordsIndex['The Effective Compliance Matrix'] + sm.eff_props[3]['compliance'] = utl.textToMatrix( + linesRead[indexCompliance + 2:indexCompliance + 2 + ems]) + except KeyError: + print('No effective compliance matrix found.') + + try: + indexConstant = keywordsIndex['The Engineering Constants'] + for line in linesRead[indexConstant + 2:indexConstant + 11]: + line = line.strip() + line = line.split('=') + label = line[0].strip().lower() + value = float(line[-1].strip()) + sm.eff_props[3]['constants'][label] = value + except KeyError: + print('No engineering constants found.') + + line = linesRead[keywordsIndex['Effective Density']] + line = line.strip().split('=') + sm.eff_props[3]['density'] = float(line[-1].strip()) + + return sm + + +def readSCOutFailure(fn_swiftcomp_in, failure_analysis): + fn_swiftcomp_in += '.fi' + lines = [] + kw_index = {} + results = [] + with open(fn_swiftcomp_in, 'r') as fin: + for i, line in enumerate(fin): + lines.append(line) + if failure_analysis == 'f': + # initial failue strength + if 'Initial Failure Strengths' in line: + kw_index['Initial Failure Strengths'] = i + elif failure_analysis == 'fe': + # initial failure envelope + pass + elif failure_analysis == 'fi': + # initial failure indices and strength ratios + pass + + if failure_analysis == 'f': + # initial failure strength + try: + index = kw_index['Initial Failure Strengths'] + for line in lines[index + 2:index + 8]: + line = line.strip().split() + results.append(map(float, line[:2])) + except KeyError: + print('No initial failure strength found.') + elif failure_analysis == 'fe': + # initial failure envelope + for line in lines: + line = line.strip().split() + results.append([int(line[0]), float(line[1]), float(line[2])]) + elif failure_analysis == 'fi': + # initial failure indices and strength ratios + for line in lines: + line = line.strip().split() + results.append([int(line[0]), float(line[1]), float(line[2])]) + + return results + + +# ==================================================================== +# +# +# WRITERS +# +# +# ==================================================================== + +def writeSCNodes(sg, fobj, sfi, sff): + if sg.sgdim == 1: + nid = sg.elements[sg.elementids1d[0]][0] + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, sg.nodes[nid]) + for eid in sg.elementids1d: + if len(sg.elements[eid]) > 2: + # node 3 + nid = sg.elements[eid][2] + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, sg.nodes[nid]) + if len(sg.elements[eid]) == 5: + # node 5 + nid = sg.elements[eid][4] + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, sg.nodes[nid]) + if len(sg.elements[eid]) > 3: + # node 4 + nid = sg.elements[eid][3] + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, sg.nodes[nid]) + # node 2 + nid = sg.elements[eid][1] + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, sg.nodes[nid]) + else: + for nid, ncoord in sg.nodes.items(): + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, ncoord) + return + + +def writeSCElements(sg, fobj, sfi): + for eid in sg.elementids1d: + elem = np.zeros(7, dtype=int) + elem[0] = eid + elem[1] = sg.elem_prop[eid] # mid/cid + for i, nid in enumerate(sg.elements[eid]): + elem[i+2] = nid + utl.writeFormatIntegers(fobj, elem) + + for eid in sg.elementids2d: + elem = np.zeros(11, dtype=int) + elem[0] = eid + elem[1] = sg.elem_prop[eid] + elem_type = 'quad' + if (len(sg.elements[eid]) == 3) or (len(sg.elements[eid]) == 6): + elem_type = 'tri' + for i, nid in enumerate(sg.elements[eid]): + if (i >= 3) and (elem_type == 'tri'): + elem[i+3] = nid + else: + elem[i+2] = nid + utl.writeFormatIntegers(fobj, elem) + + for eid in sg.elementids3d: + elem = np.zeros(22, dtype=int) + elem[0] = eid + elem[1] = sg.elem_prop[eid] + elem_type = 'hexa' + if (len(sg.elements[eid]) == 4) or (len(sg.elements[eid]) == 10): + elem_type = 'tetra' + for i, nid in enumerate(sg.elements[eid]): + if (i >= 4) and (elem_type == 'tetra'): + elem[i+3] = nid + else: + elem[i+2] = nid + utl.writeFormatIntegers(fobj, elem) + + return + + +def writeSCElementOrientations(sg, fobj, sfi, sff): + for eid, orient in sg.elem_orient.items(): + fobj.write(sfi.format(eid)) + utl.writeFormatFloats(fobj, np.hstack(orient)) + return + + +def writeSCMOCombos(sg, fobj, sfi, sff): + for cid, combo in sg.mocombos.items(): + fobj.write((sfi + sfi + sff + '\n').format(cid, combo[0], combo[1])) + return + + +def writeSCMaterials(sg, fobj, sfi, sff): + for mid, m in sg.materials.items(): + # itype = 0 + # if m.type == 'orthotropic': + # itype = 1 + # elif m.type == 'anisotropic': + # itype = 2 + mp = m.eff_props[3] + utl.writeFormatIntegers(fobj, (mid, mp['type'], 1)) + utl.writeFormatFloats(fobj, (0, mp['density'])) + if mp['type'] == 0: + mpc = mp['constants'] + utl.writeFormatFloats(fobj, [mpc['e'], mpc['nu']]) + elif mp['type'] == 1: + mpc = mp['constants'] + utl.writeFormatFloats(fobj, [mpc['e1'], mpc['e2'], mpc['e3']]) + utl.writeFormatFloats(fobj, [mpc['g12'], mpc['g13'], mpc['g23']]) + utl.writeFormatFloats(fobj, [mpc['nu12'], mpc['nu13'], mpc['nu23']]) + elif mp['type'] == 2: + for i in range(6): + for j in range(i, 6): + fobj.write(sff.format(mp['stiffness'][i][j])) + fobj.write('\n') + fobj.write('\n') + return + + +def writeSCInH(sg, fn, sfi, sff): + with open(fn, 'w') as fobj: + # Extra inputs for dimensionally reducible structures + if (sg.smdim == 1) or (sg.smdim == 2): + # model (0: classical, 1: shear refined) + fobj.write((sfi + '\n').format(sg.model)) + if sg.smdim == 1: # beam + # initial twist/curvatures + fobj.write((sff * 3 + '\n').format(0., 0., 0.)) + # oblique cross section + fobj.write((sff * 2 + '\n').format(1., 0.)) + elif sg.smdim == 2: # shell + # initial twist/curvatures + fobj.write((sff * 2 + '\n').format( + sg.initial_curvature[0], sg.initial_curvature[1] + )) + fobj.write('\n') + + # Head + fobj.write((sfi * 4 + '\n').format( + sg.physics, + sg.degen_element, + sg.trans_element, + sg.nonuniform_temperature + )) + fobj.write((sfi * 6 + '\n').format( + sg.sgdim, + len(sg.nodes), + len(sg.elements), + len(sg.materials), + sg.num_slavenodes, + len(sg.mocombos) + )) + fobj.write('\n') + + # Nodal coordinates + writeSCNodes(sg, fobj, sfi, sff) + fobj.write('\n') + + # Element connectivities + writeSCElements(sg, fobj, sfi) + fobj.write('\n') + + if sg.trans_element != 0: + writeSCElementOrientations(sg, fobj, sfi, sff) + fobj.write('\n') + + # Material-orientation combinations + if len(sg.mocombos) > 0: + writeSCMOCombos(sg, fobj, sfi, sff) + fobj.write('\n') + + # Material + writeSCMaterials(sg, fobj, sfi, sff) + fobj.write('\n') + + # Omega + fobj.write((sff + '\n').format(sg.omega)) + fobj.write('\n') + + fobj.write('\n') + + +def writeSCInD(sg, fn, sfi, sff): + with open(fn, 'w') as fobj: + utl.writeFormatFloats(fobj, sg.global_displacements, sff[2:-1]) + utl.writeFormatFloatsMatrix(fobj, sg.global_rotations, sff[2:-1]) + fobj.write((sfi+'\n').format(sg.global_loads_type)) + utl.writeFormatFloats(fobj, sg.global_loads, sff[2:-1]) + + +def writeSCInF(sg, fn, sfi, sff): + with open(fn, 'w') as fobj: + # Addtional material properties + for i, m in sg.materials.items(): + utl.writeFormatIntegers( + fobj, + (m.strength['criterion'], len(m.strength['constants'])), + sfi[2:-1] + ) + fobj.write((sff+'\n').format(m.strength['chara_len'])) + utl.writeFormatFloats(fobj, m.strength['constants'], sff[2:-1]) + + # Load type + fobj.write((sfi+'\n').format(sg.global_loads_type)) + + # (fe) Axes + + # (fi) Loads + utl.writeFormatFloats(fobj, sg.global_loads, sff[2:-1]) + + +def writeSCIn(sg, fn, analysis='h'): + """ Write SG input files for SwiftComp + + :param analysis: Type of analysis. + - 'h': Homogenization + - 'd' or 'l': Dehomogenization + - 'f': Failure analysis + """ + + if not isinstance(sg, sgm.StructureGene): + return + + # string format + sfi = '{:8d}' + sff = '{:16.6E}' + + if 'h' in analysis: + writeSCInH(sg, fn, sfi, sff) + elif ('d' in analysis) or ('l' in analysis): + writeSCInD(sg, fn, sfi, sff) + + return fn diff --git a/io/iosc.pyc b/io/iosc.pyc new file mode 100644 index 0000000..d0d4a6c Binary files /dev/null and b/io/iosc.pyc differ diff --git a/io/iovabs.py b/io/iovabs.py new file mode 100644 index 0000000..36bd573 --- /dev/null +++ b/io/iovabs.py @@ -0,0 +1,482 @@ +import os +import sys +import traceback +import numpy as np +import msgpi.io.utils as utl +import msgpi.sg as sgm + + +# ==================================================================== +# +# +# READERS +# +# +# ==================================================================== + +def readVABSIn(fn_vabs_in): + """ Read data from the VABS input file + + :param fn_vabs_in: File name of the VABS input file + :type fn_vabs_in: string + """ + + try: + fn_base, fn_extn = os.path.splitext(fn_vabs_in) + name = os.path.basename(fn_base) + sg = sgm.StructureGene(name, 2, 1) + + flag_curve = 0 + flag_oblique = 0 + + count = 0 + count_node = 1 + count_element = 1 + count_layertype = 1 + count_material = 1 + line_material = 0 + head = 3 # at least 3 lines in the head (flag) part + with open(fn_vabs_in, 'r') as fin: + for li, line in enumerate(fin): + line = line.strip() + if line == '\n' or line == '': + continue + + count += 1 + line = line.split() + if count == 1: + # format_flag nlayer + line = list(map(int, line)) + flag_format = line[0] + num_layertypes = line[1] + continue + elif count == 2: + # timoshenko_flag recover_flag thermal_flag + line = list(map(int, line)) + sg.model = line[0] + sg.analysis = line[1] + sg.physics = line[2] + if sg.physics > 0: + sg.physics == 1 + continue + elif count == 3: + # curve_flag oblique_flag trapeze_flag vlasov_flag + line = list(map(int, line)) + flag_curve = line[0] + flag_oblique = line[1] + if line[2] == 1: + sg.model = 3 # trapeze effect + if line[3] == 1: + sg.model = 2 # Vlasov model + if flag_curve == 1: + head += 1 # extra line in the head + if flag_oblique == 1: + head += 1 # extra line in the head + continue + elif flag_curve != 0 and count <= head: + pass + elif flag_oblique != 0 and count <= head: + pass + elif count == (head + 1): + # nnode nelem nmate + line = list(map(int, line)) + num_nodes = line[0] + num_elements = line[1] + num_materials = line[2] + continue + elif count_node <= num_nodes: + # Read nodal coordinates + sg.nodes[int(line[0])] = [ + float(line[1]), float(line[2])] + count_node += 1 + continue + elif count_element <= num_elements: + # Read element connectivities + eid = int(line[0]) + sg.elementids.append(eid) + temp = [int(i) for i in line[1:] if i != '0'] + sg.elements[eid] = temp + sg.elementids2d.append(eid) + count_element += 1 + continue + elif count_element <= (2 * num_elements): + # Read element theta1 + if flag_format == 1: + # new format + eid = int(line[0]) + lyrtp = int(line[1]) + theta1 = float(line[2]) + else: + # old format + eid = int(line[0]) + mid = int(line[1]) + theta3 = float(line[2]) + theta1 = float(line[3]) + lyrtp = 0 + for k, v in sg.mocombos.items(): + if (mid == v[0]) and (theta3 == v[1]): + lyrtp = k + break + if lyrtp == 0: + lyrtp = len(sg.mocombos) + 1 + sg.mocombos[lyrtp] = [mid, theta3] + + if lyrtp not in sg.prop_elem.keys(): + sg.prop_elem[lyrtp] = [] + sg.prop_elem[lyrtp].append(eid) + sg.elem_prop[eid] = lyrtp + # self.elem_angle[eid] = theta1 + sg.elem_orient[eid] = [ + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 0.0] + ] + sg.elem_orient[eid][1][1] = np.cos(theta1) + sg.elem_orient[eid][1][2] = np.sin(theta1) + count_element += 1 + continue + elif (flag_format == 1) and (count_layertype <= num_layertypes): + # Read layer types + ltid = int(line[0]) + mid = int(line[1]) + theta3 = float(line[2]) + sg.mocombos[ltid] = [mid, theta3] + count_layertype += 1 + continue + elif count_material <= num_materials: + # Read materials + if line_material == 0: + mid = int(line[0].split(',')[0]) + mtype = int(line[1]) + sg.materials[mid] = sgm.MaterialSection() + sg.materials[mid].eff_props[3]['type'] = mtype + line_material += 1 + continue + + if mtype == 0: + # isotropic + if line_material == 1: + e = float(line[0]) + nu = float(line[1]) + sg.materials[mid].eff_props[3]['constants']['e'] = e + sg.materials[mid].eff_props[3]['constants']['nu'] = nu + line_material = 99 + continue + elif mtype == 1: + # orthotropic + if line_material == 1: + # pe = [] + e = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['e1'] = e[0] + sg.materials[mid].eff_props[3]['constants']['e2'] = e[1] + sg.materials[mid].eff_props[3]['constants']['e3'] = e[2] + # for ei in e: + # pe.append(ei) + # sg.materials[mid].property_elastic += e + line_material += 1 + continue + elif line_material == 2: + g = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['g12'] = g[0] + sg.materials[mid].eff_props[3]['constants']['g13'] = g[1] + sg.materials[mid].eff_props[3]['constants']['g23'] = g[2] + # for gi in g: + # pe.append(gi) + # self.materials[mid].property_elastic += g + line_material += 1 + continue + elif line_material == 3: + nu = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['nu12'] = nu[0] + sg.materials[mid].eff_props[3]['constants']['nu13'] = nu[1] + sg.materials[mid].eff_props[3]['constants']['nu23'] = nu[2] + # for nui in nu: + # pe.append(nui) + # sg.materials[mid].property_elastic.append(pe) + line_material = 99 + continue + elif mtype == 2: + # anisotropic + continue + + if line_material == 99: + dens = float(line[0]) + sg.materials[mid].eff_props[3]['density'] = dens + # self.mate_propt[mid].append(dens) + line_material = 0 + count_material += 1 # next material + continue + continue + except: + # e = sys.exc_info()[0] + e = traceback.format_exc() + print(e) + + return sg + + +def readVABSOutHomo(fn, scrnout=True): + """Read VABS homogenization results + + :param fn: VABS output file name (e.g. example.sg.k) + :type fn: string + """ + sm = sgm.MaterialSection(1) + + linesRead = [] + keywordsIndex = {} + + with open(fn, 'r') as fin: + for lineNumber, line in enumerate(fin): + linesRead.append(line) + if 'The Mass Center of the Cross Section' in line: + keywordsIndex['mc'] = lineNumber + if 'The Geometric Center of the Cross Section' in line: + keywordsIndex['gc'] = lineNumber + if 'The Neutral Axes (or Tension Center) of the Cross Section' in line: + keywordsIndex['tc'] = lineNumber + if 'The Generalized Shear Center of the Cross Section in the User Coordinate System' in line: + keywordsIndex['sc'] = lineNumber + if 'Mass Per Unit Span' in line: + sm.eff_props[1]['density'] = float(line.split('=')[-1].strip()) + if 'Classical Stiffness Matrix' in line: + keywordsIndex['csm'] = lineNumber + if 'Classical Flexibility Matrix' in line: + keywordsIndex['cfm'] = lineNumber + if 'Timoshenko Stiffness Matrix' in line: + keywordsIndex['tsm'] = lineNumber + if 'Timoshenko Flexibility Matrix' in line: + keywordsIndex['tfm'] = lineNumber + if 'Mass Matrix' in line: + if 'Mass Center' in line: + keywordsIndex['massatcenter'] = lineNumber + else: + keywordsIndex['mass'] = lineNumber + + ln = keywordsIndex['mc'] + mc = np.array([ + 0.0, + float(linesRead[ln + 3].split()[-1]), + float(linesRead[ln + 4].split()[-1]) + ]) + + ln = keywordsIndex['gc'] + gc = np.array([ + 0.0, + float(linesRead[ln + 3].split()[-1]), + float(linesRead[ln + 4].split()[-1]) + ]) + + ln = keywordsIndex['tc'] + tc = np.array([ + 0.0, + float(linesRead[ln + 3].split()[-1]), + float(linesRead[ln + 4].split()[-1]) + ]) + + if 'sc' in keywordsIndex.keys(): + ln = keywordsIndex['sc'] + sm.eff_props[1]['shearcenter'] = [ + float(linesRead[ln + 3].split()[-1]), float(linesRead[ln + 4].split()[-1]) + ] + + ln = keywordsIndex['mass'] + sm.eff_props[1]['mass']['origin'] = utl.textToMatrix(linesRead[ln + 3:ln + 9]) + try: + ln = keywordsIndex['massatcenter'] + sm.eff_props[1]['mass']['masscenter'] = utl.textToMatrix(linesRead[ln + 3:ln + 9]) + except KeyError: + sm.eff_props[1]['mass']['masscenter'] = sm.eff_props[1]['mass']['origin'] + + ln = keywordsIndex['csm'] + sm.eff_props[1]['stiffness']['classical'] = utl.textToMatrix(linesRead[ln + 3:ln + 7]) + ln = keywordsIndex['cfm'] + sm.eff_props[1]['compliance']['classical'] = utl.textToMatrix(linesRead[ln + 3:ln + 7]) + + try: + ln = keywordsIndex['tsm'] + sm.eff_props[1]['stiffness']['refined'] = utl.textToMatrix(linesRead[ln + 3:ln + 9]) + except KeyError: + if scrnout: + print('No timoshenko stiffness matrix found.') + else: + pass + try: + ln = keywordsIndex['tfm'] + sm.eff_props[1]['compliance']['refined'] = utl.textToMatrix(linesRead[ln + 3:ln + 9]) + except KeyError: + if scrnout: + print('No timoshenko flexibility matrix found.') + else: + pass + + return sm + + +# ==================================================================== +# +# +# WRITERS +# +# +# ==================================================================== + +def writeVABSNodes(sg, fobj, sfi, sff): + for nid, ncoord in sg.nodes.items(): + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, ncoord) + return + + +def writeVABSElements(sg, fobj, sfi): + for eid in sg.elementids2d: + elem = np.zeros(10, dtype=int) + elem[0] = eid + elem_type = 'quad' + if (len(sg.elements[eid]) == 3) or (len(sg.elements[eid]) == 6): + elem_type = 'tri' + for i, nid in enumerate(sg.elements[eid]): + if (i >= 3) and (elem_type == 'tri'): + elem[i+2] = nid + else: + elem[i+1] = nid + utl.writeFormatIntegers(fobj, elem) + + return + + +def writeVABSElementOrientations(sg, fobj, sfi, sff, fmt): + for eid, orient in sg.elem_orient.items(): + fobj.write(sfi.format(eid)) + fobj.write(sfi.format(sg.elem_prop[eid])) + theta1 = np.arctan2(orient[1][2], orient[1][1]) + fobj.write(sff.format(theta1)) + fobj.write('\n') + # utl.writeFormatFloats(fobj, np.hstack(orient)) + return + + +def writeVABSMOCombos(sg, fobj, sfi, sff): + for cid, combo in sg.mocombos.items(): + fobj.write((sfi + sfi + sff + '\n').format(cid, combo[0], combo[1])) + return + + +def writeVABSMaterials(sg, fobj, sfi, sff): + for mid, m in sg.materials.items(): + # itype = 0 + # if (m.type == 'orthotropic') or (m.type == 1): + # itype = 1 + # elif (m.type == 'anisotropic') or (m.type == 2): + # itype = 2 + mp = m.eff_props[3] + utl.writeFormatIntegers(fobj, (mid, mp['type'])) + + if mp['type'] == 0: + mpc = mp['constants'] + utl.writeFormatFloats(fobj, [mpc['e'], mpc['nu']]) + elif mp['type'] == 1: + mpc = mp['constants'] + utl.writeFormatFloats(fobj, [mpc['e1'], mpc['e2'], mpc['e3']]) + utl.writeFormatFloats(fobj, [mpc['g12'], mpc['g13'], mpc['g23']]) + utl.writeFormatFloats(fobj, [mpc['nu12'], mpc['nu13'], mpc['nu23']]) + elif mp['type'] == 2: + for i in range(6): + for j in range(i, 6): + fobj.write(sff.format(mp['stiffness'][i][j])) + fobj.write('\n') + + utl.writeFormatFloats(fobj, (0, mp['density'])) + fobj.write('\n') + return + + +def writeVABSMacroData(sg, fobj, sff): + utl.writeFormatFloats(fobj, sg.macro_displacements) + fobj.write('\n') + utl.writeFormatFloats(fobj, sg.macro_rotations[0]) + utl.writeFormatFloats(fobj, sg.macro_rotations[1]) + utl.writeFormatFloats(fobj, sg.macro_rotations[2]) + fobj.write('\n') + if sg.model == 0: + utl.writeFormatFloats(fobj, sg.macro_loads) + else: + utl.writeFormatFloats(fobj, [sg.macro_loads[i] for i in [0, 3, 4, 5]]) + utl.writeFormatFloats(fobj, [sg.macro_loads[i] for i in [1, 2]]) + fobj.write('\n') + utl.writeFormatFloats(fobj, sg.macro_loads_dist[0]) + utl.writeFormatFloats(fobj, sg.macro_loads_dist[1]) + utl.writeFormatFloats(fobj, sg.macro_loads_dist[2]) + utl.writeFormatFloats(fobj, sg.macro_loads_dist[3]) + return + + +def writeVABSIn(sg, fn_vabs_in='', fmt=1): + if not isinstance(sg, sgm.StructureGene): + return + + # string format + sfi = '{:8d}' + sff = '{:16.6E}' + + with open(fn_vabs_in, 'w') as fout: + utl.writeFormatIntegers(fout, [fmt, len(sg.mocombos)]) + model = 0 + trapeze = 0 + vlasov = 0 + if sg.model == 1: + model = 1 + elif sg.model == 2: + model = 1 + vlasov = 1 + elif sg.model == 3: + trapeze = 1 + physics = 0 + if sg.physics == 1: + physics = 3 + utl.writeFormatIntegers( + fout, [model, sg.analysis, physics]) + line = [0, 0, trapeze, vlasov] + if (sg.initial_twist != 0.0) or (sg.initial_curvature[0] != 0.0) or (sg.initial_curvature[1] != 0.0): + line[0] = 1 + if (sg.oblique[0] != 1.0) or (sg.oblique[1] != 0.0): + line[1] = 1 + utl.writeFormatIntegers(fout, line) + if line[0] == 1: + utl.writeFormatFloats( + fout, + [sg.initial_twist, sg.initial_curvature[0], sg.initial_curvature[1]] + ) + if line[1] == 1: + utl.writeFormatFloats(fout, sg.oblique) + utl.writeFormatIntegers( + fout, [len(sg.nodes), len(sg.elements), len(sg.materials)]) + fout.write('\n') + + # Nodes + writeVABSNodes(sg, fout, sfi, sff) + fout.write('\n') + + # Elements + writeVABSElements(sg, fout, sfi) + fout.write('\n') + + # Element, layer type, theta 1 + writeVABSElementOrientations(sg, fout, sfi, sff, fmt) + fout.write('\n') + + if fmt == 1: + # Layer types + writeVABSMOCombos(sg, fout, sfi, sff) + fout.write('\n') + + # Materials + writeVABSMaterials(sg, fout, sfi, sff) + fout.write('\n') + + # Recover + if sg.analysis == 1: + writeVABSMacroData(sg, fout, sff) + fout.write('\n') + + return fn_vabs_in diff --git a/io/iovabs.pyc b/io/iovabs.pyc new file mode 100644 index 0000000..e7a36a9 Binary files /dev/null and b/io/iovabs.pyc differ diff --git a/io/utils.py b/io/utils.py new file mode 100644 index 0000000..aa948ad --- /dev/null +++ b/io/utils.py @@ -0,0 +1,88 @@ +import numpy as np + + +def writeFormatIntegers(file, numbers, fmt='8d', newline=True): + ''' Write a list of integers into a file + + :param file: The file object for writing + :type file: file + + :param numbers: The list of numbers that is going to be written + :type numbers: list + + :param fmt: The desired format for each number + :type fmt: string + + :param newline: If append the character ``\\n`` after writting all numbers or not + :type newline: bool + ''' + fmt = '{0:' + fmt + '}' + for i in numbers: + file.write(fmt.format(i)) + if newline: + file.write('\n') + return + + +def writeFormatFloats(file, numbers, fmt='16.6E', newline=True): + ''' Write a list of floats into a file + + :param file: The file object for writing + :type file: file + + :param numbers: The list of numbers that is going to be written + :type numbers: list + + :param fmt: The desired format for each number + :type fmt: string + + :param newline: If append the character ``\\n`` after writting all numbers or not + :type newline: bool + ''' + fmt = '{0:' + fmt + '}' + for f in numbers: + file.write(fmt.format(f)) + if newline: + file.write('\n') + return + + +def writeFormatFloatsMatrix(fobj, matrix, fmt='16.6E'): + for row in matrix: + writeFormatFloats(fobj, row, fmt) + return + + +def writeFormatIntegersMatrix(fobj, matrix, fmt='8d'): + for row in matrix: + writeFormatIntegers(fobj, row, fmt) + return + + +def textToMatrix(textList): + ''' Convert the text of a block of numbers into a matrix + + :param textList: The block of text representing a matrix + :type textList: list + + :return: A matrix of numbers + :rtype: numpy.array + + :Example: + + >>> lines = ['1 2 3', '4 5 6', '7 8 9'] + >>> utilities.textToMatrix(lines) + array([[1., 2., 3.], + [4., 5., 6.], + [7., 8., 9.]]) + ''' + matrix = [] + for line in textList: + line = line.strip() + line = line.split() + row = list(map(float, line)) + # for j in line: + # row.append(float(j)) + matrix.append(row) + + return np.array(matrix) diff --git a/io/utils.pyc b/io/utils.pyc new file mode 100644 index 0000000..5dffa6b Binary files /dev/null and b/io/utils.pyc differ diff --git a/ms/__init__.py b/ms/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/ms/__pycache__/__init__.cpython-37.pyc b/ms/__pycache__/__init__.cpython-37.pyc new file mode 100644 index 0000000..945774d Binary files /dev/null and b/ms/__pycache__/__init__.cpython-37.pyc differ diff --git a/ms/__pycache__/analysis.cpython-37.pyc b/ms/__pycache__/analysis.cpython-37.pyc new file mode 100644 index 0000000..86e02ae Binary files /dev/null and b/ms/__pycache__/analysis.cpython-37.pyc differ diff --git a/ms/__pycache__/beam.cpython-37.pyc b/ms/__pycache__/beam.cpython-37.pyc new file mode 100644 index 0000000..750cc90 Binary files /dev/null and b/ms/__pycache__/beam.cpython-37.pyc differ diff --git a/ms/__pycache__/iogebt.cpython-37.pyc b/ms/__pycache__/iogebt.cpython-37.pyc new file mode 100644 index 0000000..9bcb3a3 Binary files /dev/null and b/ms/__pycache__/iogebt.cpython-37.pyc differ diff --git a/ms/__pycache__/prebeam.cpython-37.pyc b/ms/__pycache__/prebeam.cpython-37.pyc new file mode 100644 index 0000000..fd21327 Binary files /dev/null and b/ms/__pycache__/prebeam.cpython-37.pyc differ diff --git a/ms/analysis.py b/ms/analysis.py new file mode 100644 index 0000000..e8887e0 --- /dev/null +++ b/ms/analysis.py @@ -0,0 +1,105 @@ +import os +import sys +import datetime as dt +import subprocess as sbp +import numpy as np +import msgpi.ms.prebeam as prb +import msgpi.ms.iogebt as miogebt +import msgpi.io.iovabs as miovabs + + +def solveGEBT(beam_xml, scrnout=True): + """ Carry out a global beam analysis using GEBT + """ + + # Preprocess + beam = prb.preBeam(beam_xml) + + for sid, mso in beam.sections.items(): + section = miovabs.readVABSOutHomo(mso.name+'.sg.k') + beam.sections[sid] = section + + fn_gebt_in = miogebt.writeGEBTIn(beam, beam.name+'.dat') + + # Solve + fn_gebt_out = runGEBT(fn_gebt_in, scrnout) + + # Parse results + results = miogebt.readGEBTOut(fn_gebt_out, beam) + + return results + + +def runGEBT(fn_input, scrnout=True): + cmd = ['gebt', fn_input] + cmd = ' '.join(cmd) + + try: + sys.stdout.write(' - running gebt...\n') + sys.stdout.flush() + if scrnout: + sbp.call(cmd) + else: + FNULL = open(os.devnull, 'w') + sbp.call(cmd, stdout=FNULL, stderr=sbp.STDOUT) + except: + # sys.stdout.write('failed\n') + # sys.stdout.flush() + return + + return fn_input + '.out' + + +def solvePLECS(length, compliance, x1, f1=0, f2=0, f3=0, m1=0, m2=0, m3=0): + """ Solve the static problem of a prismatic, linearly elastic, + cantilever beam. Find the three displacements and three rotations + of a point x1 given loads f1, f2, f3, m1, m2, m3 applied at the + tip. + + Equations (5.59) and (5.61) from the book + Nonlinear Composite Beam Theory by D. H. Hodges + are used. + + :param length: Total length of the beam + :type length: float + + :param compliance: The 6x6 compliance matrix of the beam cross-section + :type compliance: list of list of float + + :param x1: The location where the result is wanted. + :type x1: float + """ + + F = np.array([[f1, f2, f3, m1, m2, m3]]).T + e1tilde = np.array([ + [0.0, 0.0, 0.0], + [0.0, 0.0, -1.0], + [0.0, 1.0, 0.0] + ]) + + Fx = np.array([[f1, f2, f3]]).T + Mx = np.array([[m1, m2, m3]]).T + (length - x1) * np.dot(e1tilde, F[:3, :]) + + R = compliance[:3, :3] + Z = compliance[:3, 3:] + T = compliance[3:, 3:] + + K = np.zeros((6, 6)) + + # Equation (5.61) + K[:3, :3] = ( + x1 * R + + (x1 * length - x1**2 / 2.0) * np.dot(Z, e1tilde) - + x1**2 / 2.0 * np.dot(e1tilde, Z.T) - + (x1**2 * length / 2.0 - x1**3 / 6.0) * np.dot(e1tilde, np.dot(T, e1tilde)) + ) + K[:3, 3:] = x1 * Z - x1**2 / 2.0 * np.dot(e1tilde, T) + + # Equation (5.59) + K[3:, :3] = x1 * Z.T + (x1 * length - x1**2 / 2.0) * np.dot(T, e1tilde) + K[3:, 3:] = x1 * T + + ur = np.dot(K, F) + result = np.vstack((ur, Fx, Mx)) + + return result diff --git a/ms/beam.py b/ms/beam.py new file mode 100644 index 0000000..80ed51f --- /dev/null +++ b/ms/beam.py @@ -0,0 +1,98 @@ +import numpy as np +import msgpi.sg as sgi + + +class BeamSegment(object): + """ Class for a beam segment. + """ + + def __init__(self): + self.points = [] #: Point labels [start, stop] + self.coords = [] #: Coordinates [[x1, x2, x3], [x1, x2, x3]] + self.css = [] #: Cross-section labels [start, stop] + + self.rotate_a1 = 0.0 + self.twist = 0.0 #: Twist + # self.dihedral = 0. #: Dihedral + # self.sweep = 0. #: Sweep + self.local_frame_id = 0 #: Local frame + self.frame_id = 0 + self.curv_id = 0 + + self.num_divisions = 1 #: Number of divisions for meshing + + def calcLengthSq(self): + return ((np.array(self.coords[0]) - np.array(self.coords[1]))**2).sum() + + +class Beam(object): + """ A slender beam-like structure + + This class is mainly for the GEBT code. + """ + + def __init__(self): + self.name = '' #: Name of the beam + + # Analysis + self.analysis_type = 0 #: Analysis type (GEBT) + self.max_iteration = 1 #: Max iteration + self.num_steps = 1 #: Number of analysis steps + self.num_eigens = 0 #: Number of eigen analysis resutls + + # Design + self.angular_velocity = [] #: Angular velocity of the rotating beam + self.linear_velocity = [] #: Linear velocity of the first key point + + #: Points + #: {ptid: [x1, x2, x3], ...} + self.points = {} + + #: Beam segments + #: {bsid: BeamSegment object, ...} + self.segments = {} + + # self.divisions = [] #: Divisions of members + self.pconditions = [] #: Point conditions (B.C. and loads) + self.mconditions = [] #: Member conditions (B.C. and loads) + + #: Effective properties of cross-sections + #: {sid: MaterialSection object, ...} + self.sections = {} + + self.frames = {} #: Local self.frames + self.distrloads = [] + self.timefunctions = [] + self.initcurvatures = [] + + # Results + # self.results = None #: Results of GEBT analysis + + def echo(self): + print('') + print('Key Point Coordinates') + print('='*40) + for pid, coords in self.points.items(): + print('point {pid:4d}: {xyz[0]:16.6e} {xyz[1]:16.6e} {xyz[2]:16.6e}'.format(pid=pid, xyz=coords)) + print('') + print('Member Definition') + print('='*40) + for mid, bs in self.segments.items(): + print( + 'member {mid:4d}: {pids[0]:4d} {pids[1]:4d} {cs[0]:4d} {cs[1]:4d} {frm:4d} {ndiv:4d} {curv:4d}'.format( + mid=mid, pids=bs.points, cs=bs.css, frm=bs.frame_id, ndiv=bs.num_divisions, curv=bs.curv_id + ) + ) + print('') + + def findPtCoordByName(self, name): + for i, c in self.points.items(): + if i == name: + return c + return None + + def findSectionByName(self, name): + for i, s in self.sections.items(): + if s.name == name: + return i + return 0 diff --git a/ms/iogebt.py b/ms/iogebt.py new file mode 100644 index 0000000..0d7103d --- /dev/null +++ b/ms/iogebt.py @@ -0,0 +1,419 @@ +import numpy as np +import msgpi.ms.beam as msb +import msgpi.sg as sgi +import msgpi.io.utils as utl + + +# ==================================================================== +# +# +# READERS +# +# +# ==================================================================== + +def readGEBTIn(fn_gebt_in): + beam = msb.Beam() + + # results = {} + lines = [] + + with open(fn_gebt_in, 'r') as fin: + for ln, line in enumerate(fin): + line = line.strip() + if line == '': + continue + else: + lines.append(line.split('#')[0].strip()) + + li = 0 + + line1 = list(map(int, lines[li].split())) + beam.analysis_type = line1[0] + beam.max_iteration = line1[1] + beam.num_steps = line1[2] + if beam.analysis_type == 0: + li += 1 + else: + li += 5 + if beam.analysis_type == 3: + beam.num_eigens = int(lines[li]) + li += 1 + + line2 = list(map(int, lines[li].split())) + nkp = line2[0] + nmemb = line2[1] + # results['ncondpt'] = line2[2] + # results['nmate'] = line2[3] + # results['nframe'] = line2[4] + # results['ncondmb'] = line2[5] + # results['ndistr'] = line2[6] + # results['ntimefun'] = line2[7] + # results['ncurv'] = line2[8] + li += 1 + + # Read keypoints + for i in range(li, li + nkp): + p = lines[i].strip().split()[:4] + beam.points[int(p[0])] = list(map(float, p[1:])) + li += nkp + + # Read members + for i in range(li, li + nmemb): + data = list(map(int, lines[i].strip().split()[:8])) + bs = msb.BeamSegment() + bs.points = data[1:3] + bs.css = data[3:5] + bs.frame_id = data[5] + bs.num_divisions = data[6] + bs.curv_id = data[7] + beam.segments[data[0]] = bs + li += nmemb + + return beam + + + + + + + +# -------------------------------------------------------------------- + +def readGEBTOutNode(lines): + out = [] + for l in lines: + out += list(map(float, l.strip().split())) + return out + + + + + + + + + +def readGEBTOutStatic(fn_gebt_out, beam): + flag_analysis = beam.analysis_type + nstep = beam.num_steps + nkp = len(beam.points) + nmemb = len(beam.segments) + + results = [] # for all steps + points_results = [] + members_results = [] + + # lines = [] + + with open(fn_gebt_out, 'r') as fin: + all_lines = fin.readlines() + + i = 0 + sid = 0 # at least one step + mid = 0 + while i < len(all_lines): + line = all_lines[i].strip() + if (line == '') or ('===' in line) or ('---' in line): + i += 1 + continue + elif ('Step' in line) or ((nstep == 1) and (sid == 0)): + sid += 1 # step number + results_step = {'point': [], 'member': []} # point results, member results + mid = 0 + i += 1 + continue + elif 'Point' in line: + # Read point result + out = readGEBTOutNode(all_lines[i + 2:i + 5]) + i += 5 + results_step['point'].append(out) + # points_results.append(out) + continue + elif 'Member' in line: + # Read member result + mid += 1 + nelem = beam.segments[mid].num_divisions # number of divisions + i += 2 + member_out = [] + for j in range(nelem): + if flag_analysis == 0: + # Static + out = readGEBTOutNode(all_lines[i:i + 3]) + i += 4 + else: + # Dynamic + out = readGEBTOutNode(all_lines[i:i + 4]) + i += 5 + member_out.append(out) + results_step['member'].append(member_out) + # members_results.append(member_out) + if mid == nmemb: + results.append(results_step) + continue + i += 1 + + return results + + + + + + + +# -------------------------------------------------------------------- + +def readGEBTOutEigen(fn_gebt_out, beam): + + # if len(gebtin) == 0: + # gebtin = readGEBTIn(gebtin_name) + # gebtout_name = gebtin_name + '.out' + + nstep = beam.num_steps + nkp = len(beam.points) + nmemb = len(beam.segments) + members_in = beam.segments + # members_in = gebtin['members'] + # print 'nmemb =', nmemb + + # evei = [p1, m11, m12, ..., m1n, p2] + # p1/m1i = [x1, x2, x3, t1, t2, t3] + + with open(fn_gebt_out, 'r') as fin: + all_lines = fin.readlines() + + + + + # Read steady state results + + results_steady = [] # for all steps + points_results = [] + members_results = [] + + end_of_steady = False + + i = 0 + sid = 0 # at least one step + mid = 0 + # while i < len(all_lines): + while not end_of_steady: + line = all_lines[i].strip() + if (line == '') or ('===' in line) or ('---' in line): + i += 1 + continue + elif ('Step' in line) or ((nstep == 1) and (sid == 0)): + sid += 1 # step number + results_step = {'point': [], 'member': []} # point results, member results + mid = 0 + i += 1 + continue + elif 'Point' in line: + # Read point result + out = readGEBTOutNode(all_lines[i + 2:i + 5]) + i += 5 + results_step['point'].append(out) + # points_results.append(out) + continue + elif 'Member' in line: + # Read member result + mid += 1 + nelem = beam.segments[mid].num_divisions # number of divisions + i += 2 + member_out = [] + for j in range(nelem): + out = readGEBTOutNode(all_lines[i:i + 4]) + i += 5 + member_out.append(out) + results_step['member'].append(member_out) + # members_results.append(member_out) + if mid == nmemb: + results_steady.append(results_step) + end_of_steady = True + continue + i += 1 + + + + + # Read eigen results + all_lines = all_lines[i:] + + eva = [] # Eigenvalues: [[eva11, eva12], [eva21, eva22], ...] + eve = [] # Eigenvectors: [eve1, eve2, ...] + + eig = 0 + i = 0 + block = '' + while i < len(all_lines): + # print i + line = all_lines[i].strip() + if line == '': + i += 1 + continue + if 'Eigenvalue' in line: + eig += 1 + # print 'eigenvalue:', eig + eva.append(list(map(float, all_lines[i + 1].split()))) + evei = [[], []] + # store all Point results for an Eigenvalue temporarily + # evei_p = [] # = [[x11, x12, x13, u11, u12, u13, t11, t12, t13], ...] + # store all Member results for an Eigenvalue temporarily + # evei_m = [] # = [[[x11, x12, x13, u11, u12, u13, t11, t12, t13], [], ...], ...] + mid = 0 # Member id + pid = 0 + i += 2 + continue + if eig > 0: + if 'Point' in line: + pid += 1 + out = readGEBTOutNode(all_lines[i + 2:i + 5]) + # pr = utl.textToMatrix(all_lines[i + 2:i + 5]) # point results + # x, u, t = pr[0], pr[1][:3], pr[1][3:] + # [x11, x12, x13, u11, u12, u13, t11, t12, t13] + # eve_p.append(pr[0] + pr[1]) + i += 5 + evei[0].append(out) + continue + elif 'Member' in line: + # Read member result + mid += 1 + nelem = beam.segments[mid].num_divisions # number of divisions + # nelem = members_in[mid - 1][-2] # number of divisions + i += 2 + evei_m = [] + for j in range(nelem): + out = readGEBTOutNode(all_lines[i:i + 4]) + i += 5 + evei_m.append(out) + evei[1].append(evei_m) + # members_results.append(member_out) + if mid == nmemb: + eve.append(evei) + continue + i += 1 + + results_eigen = [eva, eve] + + return results_steady, results_eigen + + + + + + + + + +def readGEBTOut(fn_gebt_out, beam): + flag_analysis = beam.analysis_type + if flag_analysis <= 1: + return readGEBTOutStatic(fn_gebt_out, beam) + elif flag_analysis == 3: + return readGEBTOutEigen(fn_gebt_out, beam) + + + + + + + + + + + +# ==================================================================== +# +# +# WRITERS +# +# +# ==================================================================== + +def writeGEBTIn(beam, fn_gebt_in=''): + """ Write data to the GEBT input""" + + if fn_gebt_in == '': + fn_gebt_in = beam.name + '.dat' + + with open(fn_gebt_in, 'w') as fout: + utl.writeFormatIntegers( + fout, [beam.analysis_type, beam.max_iteration, beam.num_steps], '8d' + ) + if beam.analysis_type != 0: + utl.writeFormatFloats(fout, beam.angular_velocity[0]) + utl.writeFormatIntegers(fout, beam.angular_velocity[1]) + utl.writeFormatFloats(fout, beam.linear_velocity[0]) + utl.writeFormatIntegers(fout, beam.linear_velocity[1]) + if beam.analysis_type == 3: + fout.write('{0:8d}\n'.format(beam.num_eigens)) + fout.write('\n') + + line = [ + len(beam.points), len(beam.segments), len(beam.pconditions), + len(beam.sections), len(beam.frames), len(beam.mconditions), + len(beam.distrloads), len(beam.timefunctions), len(beam.initcurvatures) + ] + utl.writeFormatIntegers(fout, line, '4d') + fout.write('\n') + + # Write the block of points + for pid, coord in beam.points.items(): + fout.write('{0:4d}'.format(pid)) + utl.writeFormatFloats(fout, coord) + fout.write('\n') + + # Write the block of members + # for i, member in enumerate(members): + for bsid, bs in beam.segments.items(): + line = [ + bsid, bs.points[0], bs.points[1], bs.css[0], bs.css[1], + bs.frame_id, bs.num_divisions, bs.curv_id + ] + utl.writeFormatIntegers(fout, line, '4d') + fout.write('\n') + + # Write the block of point conditions + for condition in beam.pconditions: + fout.write('{0:4d}\n'.format(condition['point'])) + utl.writeFormatIntegers(fout, condition['components'], '16d') + utl.writeFormatFloats(fout, condition['values']) + utl.writeFormatIntegers( + fout, np.zeros(6, dtype=np.int8), '16d') + utl.writeFormatIntegers( + fout, np.zeros(6, dtype=np.int8), '16d') + fout.write('\n') + + # Write the block of cross sectional properties + # for i, member in enumerate(members): + for sid, ms in beam.sections.items(): + fout.write('{0:4d}\n'.format(sid)) + # tcm = st['tcm'] + # for row in tcm: + # utl.writeFormatFloats(fout, row) + if len(ms.eff_props[1]['compliance']['refined']) > 0: + for row in ms.eff_props[1]['compliance']['refined']: + utl.writeFormatFloats(fout, row) + else: + for i, row in enumerate(ms.eff_props[1]['compliance']['classical']): + utl.writeFormatFloats(fout, (row[0], 0, 0, row[1], row[2], row[3])) + if i == 0: + utl.writeFormatFloats(fout, np.zeros(6)) + utl.writeFormatFloats(fout, np.zeros(6)) + fout.write('\n') + if beam.analysis_type != 0: + # mass = st['mass'] + # for row in mass: + # utl.writeFormatFloats(fout, row) + for row in ms.eff_props[1]['mass']['origin']: + utl.writeFormatFloats(fout, row) + fout.write('\n') + + # Write the block of self.frames + for i, cij in beam.frames.items(): + fout.write('{0:4d}\n'.format(i)) + for row in cij: + utl.writeFormatFloats(fout, row) + fout.write('\n') + + return fn_gebt_in diff --git a/ms/prebeam.py b/ms/prebeam.py new file mode 100644 index 0000000..8d08dd8 --- /dev/null +++ b/ms/prebeam.py @@ -0,0 +1,475 @@ +import os +import numpy as np +import xml.etree.ElementTree as et +import msgpi.ms.beam as msb +import msgpi.sg as sgi +import msgpi.cross_section as sgcs +import msgpi.utils as utl + + +def preBeam(fn_beam, mode=1, sections=[]): + """ Preprocessor of GEBT + + :param fn_beam: File name of the PreGEBT main input file (.xml) + :type fn_beam: string + + :param mode: Mode of running + 1 - create SG (cross-section) input files + 2 - run SGs + 3 - create global 1D beam input file + :type mode: int + """ + + beam = msb.Beam() + + + # ---------------------------------------------------------------- + # Parse XML parametric inputs + tree = et.parse(fn_beam) + xr_beam = tree.getroot() + # beam.name = xr_beam.get('name') + beam.name = os.path.splitext(os.path.basename(fn_beam))[0] + + # self.fn_gebt_in = self.name + '_gebt.dat' + + # Read analysis settings + xe_analysis = xr_beam.find('analysis') + beam.analysis_type = int(xe_analysis.find('type').text) + beam.max_iteration = int(xe_analysis.find('max_iteration').text) + beam.num_steps = int(xe_analysis.find('num_steps').text) + if beam.analysis_type == 3: + beam.num_eigens = int(xe_analysis.find('num_eigenvalue').text) + + + # ---------------------------------------------------------------- + # Read design + xe_design = xr_beam.find('design') + inpfmt = 1 + xa_inpfmt = xe_design.get('format') + if xa_inpfmt is not None: + inpfmt = int(xa_inpfmt) + + # Read global frame a + xe_frame_a = xe_design.find('global_frame_a') + frame_a = list(map(float, xe_frame_a.text.strip().split())) + frame_a = np.array(frame_a).reshape((3, 3)) + # print(frame_a) + + if inpfmt == 1: + # Read points + pn2l = {} # point name to label dictionary + plabel = 0 + for xe_point in xe_design.findall('point'): + plabel += 1 + pname = xe_point.get('name') + pcoord = list(map(float, xe_point.text.strip().split())) + pn2l[pname] = plabel + beam.points[plabel] = pcoord + + # Read beam segments + bsn2l = {} # beam segment name to label dictionary + bslabel = 0 + fblabel = 0 # local frame b label + for xe_sgm in xe_design.findall('segment'): + bslabel += 1 + bs = msb.BeamSegment() + + # Start station + xe_start = xe_sgm.find('start') + pn = xe_start.find('location').text.strip() + bs.points.append(pn2l[pn]) + pc = beam.findPtCoordByName(pn2l[pn]) + bs.coords.append(pc) + xe_cs = xe_start.find('cross_section') + csname = xe_cs.get('name') + cslabel = beam.findSectionByName(csname) + if cslabel == 0: + cslabel = len(beam.sections) + 1 + mso = sgi.MaterialSection(1) + mso.name = csname + for s in sections: + if s.name == csname: + mso = s + beam.sections[cslabel] = mso + if xe_start.find('rotate_a1') is not None: + bs.rotate_a1 = float(xe_start.find('rotate_a1').text.strip()) + bs.css.append(cslabel) + # print(beam.sections) + + # Stop station + xe_stop = xe_sgm.find('stop') + pn = xe_stop.find('location').text.strip() + bs.points.append(pn2l[pn]) + pc = beam.findPtCoordByName(pn2l[pn]) + bs.coords.append(pc) + xe_cs = xe_stop.find('cross_section') + csname = xe_cs.get('name') + # print(csname) + cslabel = beam.findSectionByName(csname) + # print(cslabel) + if cslabel == 0: + cslabel = len(beam.sections) + 1 + mso = sgi.MaterialSection(1) + mso.name = csname + for s in sections: + if s.name == csname: + mso = s + beam.sections[cslabel] = mso + bs.css.append(cslabel) + # print(bs.css) + + nd = int(xe_sgm.find('mesh_size').text.strip()) + bs.num_divisions = nd + + # Local frame + if xe_sgm.find('local_b') is not None: + xe_b = xe_sgm.find('local_b') + if xe_b.get('method').strip() != 'global': + p12 = np.array(list(map(float, xe_b.find('p12').text.strip().split()))) + p0 = np.array(beam.points[bs.points[0]]) + p1 = np.array(beam.points[bs.points[1]]) + b1 = p1 - p0 + b1 = b1 / np.linalg.norm(b1) + b12 = p12 - p0 + b3 = np.cross(b1, b12) + b3 = b3 / np.linalg.norm(b3) + b2 = np.cross(b3, b1) + frame_b = np.vstack((b1, b2, b3)) + # print(frame_b) + # Calculate C_ij^ab + cij = np.zeros((3, 3)) + for i in range(3): + for j in range(3): + cij[i, j] = np.dot(frame_a[i], frame_b[j]) + # print(cij) + fblabel += 1 + bs.frame_id = fblabel + beam.frames[fblabel] = cij + # for i in range(3): + # a_temp = cij[i, 0] * frame_b[0] + cij[i, 1] * frame_b[1] + cij[i, 2] * frame_b[2] + # print(a_temp) + + # Twist + twist = 0.0 + if xe_sgm.find('twist') is not None: + twist = float(xe_sgm.find('twist').text.strip()) + bs.twist = twist + + beam.segments[bslabel] = bs + + if beam.analysis_type != 0: + ws = float(xe_design.find('angular_velocity').text.strip()) + vspname = xe_design.find('linear_velocity').get('at') + vsplabel = pn2l[vspname] + vs = ws * beam.points[vsplabel][0] + ws = [0., 0., ws] + wtf = [0, 0, 0] + beam.angular_velocity = [ws, wtf] + vs = [0, vs, 0] + vtf = [0, 0, 0] + beam.linear_velocity = [vs, vtf] + + # Read point conditions + # conditions = [] + for condition in xe_design.findall('condition'): + pname = condition.find('location').text + components = list(map(int, condition.find( + 'components').text.strip().split())) + values = list(map(float, condition.find('values').text.strip().split())) + beam.pconditions.append({ + 'point': pn2l[pname], + 'components': components, + 'values': values + }) + elif inpfmt == 2: + fn_cs_base = None + fn_bsl_base = None + fn_lyp_base = None + + xe_templates = xe_design.find('templates') + if xe_templates is not None: + xe_template_cs = xe_templates.find('cross_section') + if xe_template_cs is not None: + fn_cs_base = xe_template_cs.text.strip() + + xe_template_bsl = xe_templates.find('baselines') + if xe_template_bsl is not None: + fn_bsl_base = xe_template_bsl.text.strip() + + xe_template_lyp = xe_templates.find('layups') + if xe_template_lyp is not None: + fn_lyp_base = xe_template_lyp.text.strip() + + length = float(xe_design.find('length').text.strip()) + + xe_layers = xe_design.find('layers') + + xe_domain = xe_layers.find('domain') + domain_name = xe_domain.text.strip() + + layers = [] + kps = [] + + print(' - reading layers...') + for xe_layer in xe_layers.findall('layer'): + layer = {} + layer['lamina'] = xe_layer.find('lamina').text.strip() + + layer['spanrange'] = [0, length] + xe_spanrange = xe_layer.find('spanrange') + if xe_spanrange is not None: + defby = 'normalized_ends' + xa_defby = xe_spanrange.get('define') + if xa_defby is not None: + defby = xa_defby + if defby == 'normalized_ends': + [start, stop] = list(map(float, xe_spanrange.text.strip().split())) + layer['spanrange'] = [start*length, stop*length] + + layer['angle'] = 0.0 + xe_angle = xe_layer.find('angle') + if xe_angle is not None: + func = 'constant' + xa_func = xe_angle.get('function') + if xa_func is not None: + func = xa_func + if func == 'constant': + layer['angle'] = float(xe_angle.text.strip()) + elif func == 'polynomial': + xe_order = xe_angle.find('order') + xe_coeff = xe_angle.find('coefficients') + xe_ndivs = xe_angle.find('divisions') + + porder = int(xe_order.text.strip()) + coeffs = list(map(float, xe_coeff.text.strip().split())) + ndivs = int(xe_ndivs.text.strip()) + + layer['angle_func'] = utl.Polynomial2DSP(porder, coeffs) + + layer['kps'] = np.linspace( + layer['spanrange'][0], + layer['spanrange'][1], + ndivs+1 + ) + pass + + # kps.append(layer['spanrange'][0]) + # kps.append(layer['spanrange'][1]) + kps.append(layer['kps']) + + # layer['sectionrange'] = [] + # xe_sectionrange = xe_layer.find('sectionrange') + # if xe_sectionrange is not None: + # defby = 'domain_list' + # xa_defby = xe_sectionrange.get('define') + # if xa_defby is not None: + # defby = xa_defby + # if defby == 'domain_list': + # # + + layers.append(layer) + + # for layer in layers: + # print(layer) + + kps = np.concatenate(kps).tolist() + kps = list(set(kps)) + print('kps:', kps) + + for i, kp in enumerate(kps): + xe_point = et.SubElement(xe_design, 'point') + xe_point.set('name', 'p'+str(i)) + xe_point.text = str(kp) + ' 0 0' + + # Generate beam segments + print(' - creating beam segments...') + bm_sgms = [] + layups = [] + for i, kp in enumerate(kps): + if i == 0: + # first station + bm_sgm = {} + bm_sgm['loc_start'] = 'p' + str(i) + layup = [] + for j, layer in enumerate(layers): + if kp == layer['spanrange'][0]: + layup.append(layer) + bm_sgm['layup_start'] = layup + bm_sgms.append(bm_sgm) + elif i == len(kps) - 1: + # last station + layup = [] + bm_sgms[-1]['loc_stop'] = 'p' + str(i) + for j, layer in enumerate(layers): + if kp == layer['spanrange'][1]: + layup.append(layer) + bm_sgms[-1]['layup_stop'] = layup + else: + # inner stations + layup = [] + bm_sgms[-1]['loc_stop'] = 'p' + str(i) + for j, layer in enumerate(layers): + if (kp > layer['spanrange'][0]) and (kp <= layer['spanrange'][1]): + layup.append(layer) + bm_sgms[-1]['layup_stop'] = layup + + bm_sgm = {} + bm_sgm['loc_start'] = 'p' + str(i) + layup = [] + for j, layer in enumerate(layers): + if (kp >= layer['spanrange'][0]) and (kp < layer['spanrange'][1]): + layup.append(layer) + bm_sgm['layup_start'] = layup + bm_sgms.append(bm_sgm) + + # for i, bs in enumerate(bm_sgms): + # print('') + # print('- beam segment {0}'.format(i)) + # print(' - start layup') + # for j, layer in enumerate(bs['layup_start']): + # print(j, layer) + # print(' - stop layup') + # for j, layer in enumerate(bs['layup_stop']): + # print(j, layer) + + # Write PreVABS input files + print(' - writing sg input files...') + fn_lyp_tmp = fn_lyp_base + '.xml' + fn_lyp = beam.name + '_lyps.xml' + + xt_layups = et.parse(fn_lyp_tmp) + xr_layups = xt_layups.getroot() + for i, bs in enumerate(bm_sgms): + # start layup + layup_name = '_'.join(('lyp', str(i+1), '0')) + xe_layup = et.SubElement(xr_layups, 'layup') + xe_layup.set('name', layup_name) + for j, lyr in enumerate(bs['layup_start']): + xe_layer = et.SubElement(xe_layup, 'layer') + xe_layer.set('lamina', lyr['lamina']) + + f = lyr['angle_func'] + x0, x1 = lyr['spanrange'] + y = (kps[i] - x0) / (x1 - x0) + # xe_layer.text = str(lyr['angle']) + xe_layer.text = str(f([y, 0])) + + # start cross-section + fn_cs = '_'.join((beam.name, 'cs', str(i+1), '0')) + xt_cs = et.parse(fn_cs_base+'.xml') + xr_cs = xt_cs.getroot() + xr_cs.set('name', fn_cs) + xe_include = xr_cs.find('include') + xe_include_layup = xe_include.find('layup') + xe_include_layup.text = fn_lyp[:-4] + + xe_cmp = None + for xe in xr_cs.findall('component'): + if xe.get('name') == domain_name: + xe_cmp = xe + break + if xe_cmp is None: + xe_cmp = et.SubElement(xr_cs, 'component') + xe_sgm = xe_cmp.find('segment') + if xe_sgm is None: + xe_sgm = et.SubElement(xe_cmp, 'segment') + xe_lyp = xe_sgm.find('layup') + if xe_lyp is None: + xe_lyp = et.SubElement(xe_sgm, 'layup') + xe_lyp.text = layup_name + + bs['cs_start'] = fn_cs + xt_cs.write(fn_cs+'.xml') + + # stop layup + layup_name = '_'.join(('lyp', str(i+1), '1')) + xe_layup = et.SubElement(xr_layups, 'layup') + xe_layup.set('name', layup_name) + for j, lyr in enumerate(bs['layup_stop']): + xe_layer = et.SubElement(xe_layup, 'layer') + xe_layer.set('lamina', lyr['lamina']) + + f = lyr['angle_func'] + x0, x1 = lyr['spanrange'] + y = (kps[i+1] - x0) / (x1 - x0) + # xe_layer.text = str(lyr['angle']) + xe_layer.text = str(f([y, 0])) + + # stop cross-section + fn_cs = '_'.join((beam.name, 'cs', str(i+1), '1')) + xt_cs = et.parse(fn_cs_base+'.xml') + xr_cs = xt_cs.getroot() + xr_cs.set('name', fn_cs) + xe_include = xr_cs.find('include') + xe_include_layup = xe_include.find('layup') + xe_include_layup.text = fn_lyp[:-4] + + xe_cmp = None + for xe in xr_cs.findall('component'): + if xe.get('name') == domain_name: + xe_cmp = xe + break + if xe_cmp is None: + xe_cmp = et.SubElement(xr_cs, 'component') + xe_sgm = xe_cmp.find('segment') + if xe_sgm is None: + xe_sgm = et.SubElement(xe_cmp, 'segment') + xe_lyp = xe_sgm.find('layup') + if xe_lyp is None: + xe_lyp = et.SubElement(xe_sgm, 'layup') + xe_lyp.text = layup_name + + bs['cs_stop'] = fn_cs + xt_cs.write(fn_cs+'.xml') + + # for i, bs in enumerate(bm_sgms): + # print('') + # print('- beam segment {0}'.format(i)) + # print(bs) + + xt_layups.write(fn_lyp) + + # Write beam segments + fn_beam_sgn = fn_beam[:-4] + '_sgn.xml' + print(' - writing beam segments in file ' + fn_beam_sgn) + for i, bs in enumerate(bm_sgms): + xe_bs = et.SubElement(xe_design, 'segment') + xe_start = et.SubElement(xe_bs, 'start') + xe_loc = et.SubElement(xe_start, 'location') + xe_loc.text = bs['loc_start'] + xe_cs = et.SubElement(xe_start, 'cross_section') + xe_cs.set('name', bs['cs_start']) + + xe_stop = et.SubElement(xe_bs, 'stop') + xe_loc = et.SubElement(xe_stop, 'location') + xe_loc.text = bs['loc_stop'] + xe_cs = et.SubElement(xe_stop, 'cross_section') + xe_cs.set('name', bs['cs_stop']) + + xe_design.set('format', '1') + xe_design.remove(xe_templates) + xe_design.remove(xe_layers) + tree.write(fn_beam_sgn) + + # Calculate C_ab + # self.frames = [] + frame_a = [ + np.array([1., 0., 0.]), + np.array([0., 1., 0.]), + np.array([0., 0., 1.]) + ] + # for member in members: + # for i in range(1, self.num_points): + # frame_b = self.stations[i]['frame_b'] + # Cab = utl.calcCab(frame_a, frame_b) + # self.frames.append(Cab) + # self.num_frames = len(self.frames) + + # ms = float(design.find('mesh_size').text) + # # divisions = [] + # for i in range(1, self.num_points): + # dx1 = self.stations[i]['coordinates'][0] - \ + # self.stations[i - 1]['coordinates'][0] + # self.divisions.append(int(np.ceil(dx1 / ms))) + + return beam diff --git a/presg.py b/presg.py new file mode 100644 index 0000000..47309aa --- /dev/null +++ b/presg.py @@ -0,0 +1,329 @@ +import os +import sys +import subprocess as sbp +import numpy as np +import xml.etree.ElementTree as et +import msgpi.sg as sgm +import msgpi.io.iosc as miosc +import msgpi.io.iovabs as miovabs +import msgpi.utils as utl + + +def readMaterialFromXMLElement(xem): + type_flag = { + 'isotropic': 0, + 'orthotropic': 1, + 'anisotropic': 2 + } + elastic_label = { + 0: ('e', 'nu'), + 1: ( + 'e1', 'e2', 'e3', + 'g12', 'g13', 'g23', + 'nu12', 'nu13', 'nu23' + )#, + # 'anisotropic': ( + # ('c11', 'c12', 'c13', 'c14', 'c15', 'c16'), + # ('c12', 'c22', 'c23', 'c24', 'c25', 'c26'), + # ('c13', 'c23', 'c33', 'c34', 'c35', 'c36'), + # ('c14', 'c24', 'c34', 'c44', 'c45', 'c46'), + # ('c15', 'c25', 'c35', 'c45', 'c55', 'c56'), + # ('c16', 'c26', 'c36', 'c46', 'c56', 'c66'), + # ) + } + + m = sgm.MaterialSection() + m.name = xem.get('name').strip() + + m.eff_props[3]['type'] = type_flag[xem.get('type').strip()] + m.eff_props[3]['density'] = float(xem.find('density').text.strip()) + + ep = xem.find('elastic') + if (m.eff_props[3]['type'] == 0) or (m.eff_props[3]['type'] == 1): + ep_labels = elastic_label[m.eff_props[3]['type']] + for tag in ep_labels: + m.eff_props[3]['constants'][tag] = float( + ep.find(tag).text.strip() + ) + else: + stf = np.zeros((6, 6)) + for i in range(6): + rowtag = 'row' + str(i) + row = list(map(float, ep.find(rowtag).text.strip().split())) + stf[i, i:] = row + for i in range(1, 6): + for j in range(0, i): + stf[i, j] = stf[j, i] + m.eff_props[3]['stiffness'] = stf + + if xem.find('strength') is not None: + xe_str = xem.find('strength') + if xe_str.find('criterion') is not None: + m.strength['criterion'] = int(xe_str.find('criterion').text.strip()) + m.strength['constants'] = list(map(float, xe_str.find('constants').text.strip().split())) + + return m + + +def preSG1D(xr_sg, smdim): + # debug = True + # print('running presc...') + + # ==================================================================== + # Parse XML parametric inputs + # tree = et.parse(fn_xml) + # xr_sg = sg_xml_tree.getroot() + name = xr_sg.get('name') + # sgdim = int(xr_sg.get('dim').text) + # smdim = int(xr_sg.get('model_dim').text) + + sg = sgm.StructureGene(name, 1, smdim) + + # Head + sg.analysis = int(xr_sg.find('analysis').text) + sg.model = int(xr_sg.find('model').text) # model (0: classical, 1: shear refined) + sg.degen_element = int(xr_sg.find('elem_flag').text) + sg.trans_element = int(xr_sg.find('trans_flag').text) + sg.uniform_temperature = int(xr_sg.find('temp_flag').text) + + if smdim == 2: + if xr_sg.find('init_curv') is not None: + ic = xr_sg.find('init_curv').text.strip() + sg.initial_curvature = list(map(float, ic.split())) + + # Materials and Laminae + materials = {} + d_lamina = {} + + if xr_sg.find('include') is not None: + xe_include = xr_sg.find('include') + if xe_include.find('material') is not None: + fn_material = xe_include.find('material').text.strip() + '.xml' + xt = et.parse(fn_material) + xe_materials = xt.getroot() + for xem in xe_materials.findall('material'): + m = readMaterialFromXMLElement(xem) + materials[m.name] = m + for lmn in xe_materials.findall('lamina'): + ln = lmn.get('name').strip() + lm = lmn.find('material').text.strip() + lt = float(lmn.find('thickness').text.strip()) + d_lamina[ln] = {'material': lm, 'thickness': lt} + + if xr_sg.find('materials') is not None: + xe_materials = xr_sg.find('materials') + for xem in xe_materials.findall('material'): + m = readMaterialFromXMLElement(xem) + materials[m.name] = m + for lmn in xe_materials.findall('lamina'): + ln = lmn.get('name').strip() + lm = lmn.find('material').text.strip() + lt = float(lmn.find('thickness').text.strip()) + d_lamina[ln] = {'material': lm, 'thickness': lt} + + # Layup + # List of material, thickness, angle + # [{'material': 'name', 'thickness': t, 'angle': a, 'combo': id}, {}, ...] + ld_layer = [] + l_material_used = [] + layup = xr_sg.find('layup') + method = 'll' # ll: layer list, ss: stacking sequence, pp: ply percentage + if layup.get('method') is not None: + method = layup.get('method') + if method == 'll': + # layer list + for lyr in layup.findall('layer'): + ln = lyr.get('lamina').strip() + lm = d_lamina[ln]['material'] + # if lm not in l_material_used: + # l_material_used.append(lm) + # d_material[lm]['id'] = len(l_material_used) + lt = d_lamina[ln]['thickness'] + nl = 1 + lon = lyr.text.strip().split(':') + lo = float(lon[0].strip()) + if len(lon) > 1: + nl = int(lon[1].strip()) + for ni in range(nl): + ld_layer.append( + {'material': lm, 'thickness': lt, 'angle': lo} + ) + elif method == 'ss': + # stacking sequence + ln = layup.find('lamina').text.strip() + lm = d_lamina[ln]['material'] + # l_material_used.append(lm) + # d_material[lm]['id'] = len(l_material_used) + lt = d_lamina[ln]['thickness'] + layup_code = layup.find('code').text.strip() + code = utl.parseLayupCode(layup_code) + for lo in code: + ld_layer.append({'material': lm, 'thickness': lt, 'angle': lo}) + elif method == 'pp': + # ply percentage + pass + # print code + + # Record used materials and create material-orientation combinations + mid = 0 + cid = 0 + tt = 0.0 # total thickness + for layer in ld_layer: + # Check used material-orientation combination first + cid = sg.findComboByMaterialOrientation(layer['material'], layer['angle']) + if cid == 0: + # Not found + # Check used material + mid = sg.findMaterialByName(layer['material']) + if mid == 0: + # Not found + mid = len(sg.materials) + 1 + sg.materials[mid] = materials[layer['material']] + cid = len(sg.mocombos) + 1 + sg.mocombos[cid] = [mid, layer['angle']] + sg.prop_elem[cid] = [] + layer['mocombo'] = cid + tt = tt + layer['thickness'] + + # ---------------------------------------------------------------- + # Meshing + pan = 0.0 + if smdim == 3: + ens = 2 # by default, use 2-node element for 3D structure + sg.omega = tt + elif smdim == 2: + ens = 5 # by default, use 5-node element for 2D structure + + if xr_sg.find('translate') is not None: + pan = float(xr_sg.find('translate').text.strip()) + + gms = float(xr_sg.find('global_mesh_size').text.strip()) + + if xr_sg.find('element_nodes') is not None: + ens = int(xr_sg.find('element_nodes').text.strip()) + + nply = len(ld_layer) + # tthk = thickness * nply + ht = tt / 2.0 + # print tthk + + # nodes_major = np.array([-ht + pan, ]) + # nid = 1 + nid1 = 1 + eid = 0 + yprev = -ht + pan + sg.nodes[nid1] = [yprev, ] + for lyr in ld_layer: + t = lyr['thickness'] + ne = int(round(t / gms)) # number of element for this layer + # lyr['nelem'] = ne + ns = np.linspace(yprev, yprev+t, ne+1) + for i in range(ne): + # nid = nid + ens - 1 + nid2 = nid1 + ens - 1 + sg.nodes[nid2] = [ns[i+1], ] + eid = eid + 1 + sg.elements[eid] = [nid1, nid2] + sg.elementids1d.append(eid) + sg.elem_prop[eid] = lyr['mocombo'] + sg.prop_elem[lyr['mocombo']].append(eid) + nid1 = nid2 + yprev = yprev + t + + # Change the order of each element + if ens > 2: + for eid in sg.elements.keys(): + nid1 = sg.elements[eid][0] + nid2 = sg.elements[eid][1] + nq0y3 = sg.nodes[nid1][0] + nq4y3 = sg.nodes[nid2][0] + nid3 = nid1 + 1 + sg.elements[eid].append(nid3) + if ens == 4: + pass + else: + nq2y3 = (nq0y3 + nq4y3) / 2.0 + nq1y3 = (nq0y3 + nq2y3) / 2.0 + nq3y3 = (nq2y3 + nq4y3) / 2.0 + if ens == 5: + # nid = nid + 1 + nid5 = nid3 + 1 + nid4 = nid5 + 1 + sg.nodes[nid3] = [nq1y3, ] + sg.nodes[nid5] = [nq2y3, ] + sg.nodes[nid4] = [nq3y3, ] + # sg.elements[eid].append(nid3) + sg.elements[eid].append(nid4) + sg.elements[eid].append(nid5) + else: + sg.nodes[nid3] = [nq2y3, ] + + return sg + + +def preSG(sg_xml, analysis, solver='swiftcomp', write_input=True, scrnout=True): + """ Preprocessor of a structure gene + + Parameters + ---------- + sg_xml : str + File name of SG design parameters (XML format) + analysis : str + The analysis that will be carried out after the generation of + the input file. + h - homogenization + d - dehomogenization/localization/recover + f - initial failure strength + fe - initial failure envelope + fi - initial failure indices and strength ratios + solver : str + Format of the generated input file ('vabs' or 'swiftcomp') + """ + if scrnout: + print(' - preprocessing structure gene...') + # fn_head, fn_tail = os.path.split(xml_sg) + fn_base, fn_extn = os.path.splitext(sg_xml) + fn_sg_in = fn_base + '.sg' + + xt = et.parse(sg_xml) + xr_sg = xt.getroot() + + if xr_sg.tag == 'cross_section': + smdim = 1 + cmd = ['prevabs', '-i', sg_xml] + if 'vabs' in solver: + cmd.append('-vabs') + elif 'swiftcomp' in solver: + cmd.append('-sc') + cmd.append('-' + analysis) + + cmd = ' '.join(cmd) + try: + FNULL = open(os.devnull, 'w') + sbp.call(cmd, stdout=FNULL, stderr=sbp.STDOUT) + return fn_sg_in, smdim + except: + e = sys.exc_info()[0] + print('[ERROR] ', e) + return + else: + try: + sgdim = int(xr_sg.get('sgdim')) + except TypeError: + print('[ERROR] Unable to read structure gene dimension (sgdim):', xr_sg.get('sgdim')) + return + + try: + smdim = int(xr_sg.get('smdim')) + except TypeError: + print('[ERROR] Unable to read structural model dimension (smdim):', xr_sg.get('smdim')) + return + + if sgdim == 1: + sg = preSG1D(xr_sg, smdim) + + if write_input: + miosc.writeSCIn(sg, fn_sg_in) + return fn_sg_in, smdim + else: + return sg diff --git a/presg.pyc b/presg.pyc new file mode 100644 index 0000000..581cae7 Binary files /dev/null and b/presg.pyc differ diff --git a/process.py b/process.py new file mode 100644 index 0000000..b63aea7 --- /dev/null +++ b/process.py @@ -0,0 +1,81 @@ +import sys +import datetime as dt +import msgpi.analysis as sga +import dakota.interfacing as di + +# Format string +fse = 'x [{0:%H}:{0:%M}] eval {1:6d} {2:s}' +fsi = '> [{0:%H}:{0:%M}] eval {1:6d} {2:s}' + + +def printInfo(evid, *argv): + print('> [{0:%H}:{0:%M}] eval {1:6d} {2:s}'.format( + dt.datetime.now(), evid, ' '.join(argv) + )) + + +def interface(params, results, tmp_map, res_map, sg_xml, analysis, solver, scrnout=True): + """ Create a workflow for a design problem using Dakota + + Parameters + ---------- + params : Parameters + Dakota Parameters object + results : Results + Dakota Results object + tmp_map : dict + Dictionary of file names that will be used in dprepro + {'template name': 'true name', ...} + res_map : dict + Disctionary of keywords for results of Dakota + {'dakota_response_keyword': 'sg_result_keyword', ...} + sg_xml : str + Main file name for preprocessor + analysis : str + The analysis that will be carried out after the generation of + the input file. + h - homogenization + d - dehomogenization/localization/recover + f - initial failure strength + fe - initial failure envelope + fi - initial failure indices and strength ratios + solver : str + The name of the executable ('vabs' or 'swiftcomp') + """ + evid = int(params.eval_id) + + # try: + print(' - eval {evid:d}: writing input files...'.format(evid=evid)) + for fn_t, fn_i in tmp_map.items(): + di.dprepro(template=fn_t, parameters=params, output=fn_i) + + # Solve homogenization + print(' - eval {evid:d}: running solver...'.format(evid=evid)) + sm = sga.solve(sg_xml, analysis, solver, scrnout) + + # Extract results needed + print(' - eval {evid:d}: extracting results...'.format(evid=evid)) + for k, v in res_map.items(): + tmp = sm.eff_props + v = v.strip().split('.') + for i, tag in enumerate(v): + if '#' in tag: + tag = int(tag.replace('#', '')) + if i > 0: + tag = tag - 1 + tmp = tmp[tag] + results[k].function = tmp + + # except: + # e = sys.exc_info()[0] + # print('Error:', e) + # results.fail() + # with open('results.out', 'w') as fout: + # results.write(stream=fout) + # # printInfo(evid, '>>> FAIL <<<') + # return + + # else: + # results.write() + # # printInfo(evid, '>>> FINISH <<<') + diff --git a/sg.py b/sg.py new file mode 100644 index 0000000..f0780ee --- /dev/null +++ b/sg.py @@ -0,0 +1,241 @@ +import numpy as np + + +# class Material(object): +# """ Material class +# """ + +# def __init__(self): +# self.name = '' #: Name of the material +# self.density = 0. #: Density of the material +# self.type = None #: Type of the material (isotropic, orthotropic, anisotropic) + +# # Elastic properties +# # isotropic: [[e, nu]] +# # orthotropic: [[e1, e2, e3, g12, g13, g23, nu12, nu13, nu23]] +# # aniostropic: 6x6 matrix +# self.property_elastic = [] + + +class MaterialSection(object): + """ A macroscopic structure model + + 3D constitutive + 2D plate/shell + 1D beam + """ + def __init__(self, smdim=3): + self.smdim = smdim + self.name = '' + + self.eff_props = { + 3: { + 'density': 0, + 'type': 0, # 0: isotropic, 1: orthotropic, 2: anisotropic + 'stiffness': [], + 'compliance': [], + 'constants': {} + }, + 2: { + 'density': 0, + 'stiffness': { + 'classical': [], + 'refined': [] + }, + 'compliance': { + 'classical': [], + 'refined': [] + }, + 'constants': { + 'inplane': {}, + 'flexural': {} + } + }, + 1: { + 'density': 0, + 'mass': { + 'origin': [], + 'masscenter': [] + }, + 'stiffness': { + 'classical': [], + 'refined': [] + }, + 'compliance': { + 'classical': [], + 'refined': [] + }, + 'shearcenter': [] + } + } + + self.strength = { + 'criterion': 0, + 'chara_len': 0, + 'constants': [] + } + + def __str__(self): + s = '\n' + s = s + 'Effective properties of the SG\n' + s = s + 'Structure model dimension: {0}\n'.format(self.smdim) + if self.smdim == 3: + pass + elif self.smdim == 2: + pass + elif self.smdim == 1: + pass + + return s + + def summary(self): + print('') + print('Effective properties of the SG') + print('Structure model dimension: {0}'.format(self.smdim)) + + ep = self.eff_props[self.smdim] + if self.smdim == 3: + pass + elif self.smdim == 2: + pass + elif self.smdim == 1: + stf = ep['stiffness'] + print('The Effective Stiffness Matrix') + for row in stf['classical']: + print(row) + if len(stf['refined']) > 0: + print('Generalized Timoshenko Stiffness') + for row in stf['refined']: + print(row) + + +class StructureGene(object): + """ An FE level structure gene model in the theory of MSG + """ + + def __init__(self, name, sgdim, smdim): + self.name = name #: Name of the SG + self.sgdim = sgdim #: Microscopic dimension of the SG (1, 2, 3) + #: Macroscopic structure (1: beam, 2: plate/shell, 3: solid/block) + self.smdim = smdim + + self.fn_gmsh_msh = self.name + '.msh' #: File name of the Gmsh mesh file + + # Analysis configurations + #: What analysis + #: 0 - homogenization (default) + #: 1 - dehomogenization/localization/recover + #: 2 - failure (SwiftComp only) + self.analysis = 0 + + #: Physics included in the analysis + #: 0 - elastic (default) + #: 1 - thermoelastic + #: 2 - conduction + #: 3 - piezoelectric/piezomagnetic + #: 4 - thermopiezoelectric/thermopiezomagnetic + #: 5 - piezoelectromagnetic + #: 6 - thermopiezoelectromagnetic + self.physics = 0 + + #: Macroscopic structural model + #: 0 - classical (default) + #: 1 - refined (e.g. generalized Timoshenko) + #: 2 - Vlasov model (beam only) + #: 3 - trapeze effect (beam only) + self.model = 0 + + self.trans_element = 0 #: Flag of transformation of elements + self.nonuniform_temperature = 0 #: Flag of uniform temperature + + self.initial_twist = 0.0 #: Initial twist (beam only) + self.initial_curvature = [0.0, 0.0] #: Initial curvature + self.oblique = [1.0, 0.0] #: Oblique (beam only) + + # Materials + # self.num_materials = 0 #: Number of materials + self.materials = {} #: Materials {mid: MaterialSection object, ...} + + # self.num_mocombos = 0 #: Number of material-orientation combinations + #: Material-orientation combinations + #: {cid: [mid, orientation], ...} + self.mocombos = {} + + # Discretization + # self.global_mesh_size = 0 #: Global mesh size + self.degen_element = 0 #: Flag of the type of elements (SC) + # self.num_nodes = 0 #: Number of nodes + # self.num_elements = 0 #: Number of elements + self.num_slavenodes = 0 #: Number of slave nodes + + #: Nodal coordinates + #: 3D SG: {nid: [y1, y2, y3], ...} + #: 2D SG: {nid: [y2, y3], ...} + #: 1D SG: {nid: [y3], ...} + self.nodes = {} + #: Elemental connectivities {eid: [nid1, nid2, ...], ...}, no zeros + self.elements = {} + self.elementids = [] + self.elementids1d = [] + self.elementids2d = [] + self.elementids3d = [] + + #: Assignment of property to each element + self.elem_prop = {} #: {eid: mid/cid, ...} + self.prop_elem = {} #: {mid/cid: [eid, ...], ...} + + #: Element local orientations + self.elem_orient = {} #: {eid: [[a1, a2, a3], [b1, b2, b3], [c1, c2, c3]], ...} + + self.omega = 1 + + # Dehomogenization + self.global_displacements = [] #: [u1, u2, u3] + self.global_rotations = [] #: [[C11, C12, C13], [C21, C22, C23], [C31, C32, C33]] + self.global_loads_type = 0 #: 0 - generalized stresses, 1 - generalized strains + + #: Global loads + #: 3D structures + #: generalized stresses = [s11, s22, s33, s23, s13, s12] + #: generalized strains = [e11, e22, e33, e23, e13, e12] + #: Kirchhoff-Love plate/shell model + #: generalized stresses = [N11, N22, N12, M11, M22, M12] + #: generalized strains = [e11, e22, 2e12, k11, k22, 2k12] + #: Reissner-Mindlin plate/shell model + #: generalized stresses = [N11, N22, N12, M11, M22, M12, N13, N23] + #: generalized strains = [e11, e22, 2e12, k11, k22, 2k12, g13, g23] + #: Euler-Bernoulli beam model + #: generalized stresses = [F1, M1, M2, M3] + #: generalized strains = [e11, k11, k12, k13] + #: Timoshenko beam model + #: generalized stresses = [F1, F2, F3, M1, M2, M3] + #: generalized strains = [e11, g12, g13, k11, k12, k13] + self.global_loads = [] + + self.global_loads_dist = [] #: Distributed loads for Timoshenko beam model (VABS only) + + def summary(self): + print('') + print('SUMMARY OF THE SG') + print('') + print('Structure gene dimension: {0}'.format(self.sgdim)) + print('Structure model dimension: {0}'.format(self.smdim)) + print('') + print('Number of nodes: {0}'.format(len(self.nodes))) + print('Number of elements: {0}'.format(len(self.elements))) + print('') + print('Number of materials: {0}'.format(len(self.materials))) + print('Number of material-orientation combinations: {0}'.format(len(self.mocombos))) + print('') + + def findMaterialByName(self, name): + for i, m in self.materials.items(): + if m.name == name: + return i + return 0 + + def findComboByMaterialOrientation(self, name, angle): + for i, mo in self.mocombos.items(): + if (self.materials[mo[0]].name == name) and (mo[1] == angle): + return i + return 0 diff --git a/sg.pyc b/sg.pyc new file mode 100644 index 0000000..7a81f68 Binary files /dev/null and b/sg.pyc differ diff --git a/utils.py b/utils.py new file mode 100644 index 0000000..e815fef --- /dev/null +++ b/utils.py @@ -0,0 +1,281 @@ +import os +import math +import numpy as np +import pandas as pd +import xml.etree.ElementTree as et +import msgpi.sg as sgm + + +def tilde(v): + """ + ~v_{ij} = -e_{ijk} * v_k + + :param v: 3x1 vector + :type v: list or array-like + """ + vtilde = np.zeros((3, 3)) + vtilde[0, 1] = -v[2, 0] + vtilde[0, 2] = v[1, 0] + vtilde[1, 0] = v[2, 0] + vtilde[1, 2] = -v[0, 0] + vtilde[2, 0] = -v[1, 0] + vtilde[2, 1] = v[0, 0] + + return vtilde + + +def calcRotationTensorFromParameters(rp, method=''): + """ Calculate rotation tensor from rotation parameters + + :param rp: Rotation parameters (3x1) + :type rp: list-like + + :param method: Rotation parameter type + er - Eular-Rodrigues + wm - Wiener-Milenkovic + :type method: string + """ + C = np.zeros((3, 3)) + if method == 'er': + # Eular-Rodrigues + t = np.dot(rp.T, rp)[0, 0] + nmr = (1 - t / 4.0) * np.eye(3) + np.dot(rp, rp.T) / 2.0 - tilde(rp) + dnm = 1 + t / 4.0 + C = nmr / dnm + elif method == 'wm': + # Wiener-Milenkovic + t = np.dot(rp.T, rp)[0, 0] + c0 = 2 - t / 8.0 + nmr = (c0**2 - t) * np.eye(3) - 2.0 * c0 * tilde(rp) + 2.0 * np.dot(rp, rp.T) + dnm = (4 - c0)**2 + C = nmr / dnm + + return C + + +def listToString(flist, delimiter='', fmt=''): + sfmt = '{0:' + fmt + '}' + s = '' + for i, n in enumerate(flist): + if i > 0: + s = s + delimiter + ' ' + s += sfmt.format(n) + return s + + +def angleToCosine2D(angle): + angle = math.radians(angle) + c = [ + [math.cos(angle), -math.sin(angle)], + [math.sin(angle), math.cos(angle)] + ] + return c + + +def calcBasicRotation3D(angle, axis=1): + c = np.array([ + [1., 0., 0.], + [0., np.cos(angle), -np.sin(angle)], + [0., np.sin(angle), np.cos(angle)] + ]) + if axis == 2: + p = np.array([ + [0., 1., 0.], + [0., 0., 1.], + [1., 0., 0.] + ]) + c = np.dot(p.T, np.dot(c, p)) + elif axis == 3: + p = np.array([ + [0., 0., 1.], + [1., 0., 0.], + [0., 1., 0.] + ]) + c = np.dot(p.T, np.dot(c, p)) + + return c + + +def calcGeneralRotation3D(angles, order=[1, 2, 3]): + ''' Calculate the general rotation matrix + + :param angles: Three rotation angles in radians + :type angles: list + + :param order: The order of axis of the rotation operation + :type order: list + + :return: The general rotation matrix + :rtype: numpy.array + ''' + c = np.eye(3) + for angle, axis in zip(angles, order): + ci = calcBasicRotation3D(angle, axis) + c = np.dot(ci, c) + return c + + +def rotateVectorByAngle2D(v2d, angle): + c = angleToCosine2D(angle) + vp = [ + c[0][0] * v2d[0] + c[0][1] * v2d[1], + c[1][0] * v2d[0] + c[1][1] * v2d[1] + ] + return vp + + +def calcCab(a, b): + '''Calculate the direction cosine matrix between frame a and b + + :math:`C_{ij} = a_i\ \cdot\ b_j` + + :param a: list of three a basis (a_1, a_2, a_3) + :type a: list + + :param b: list of three b basis (b_1, b_2, b_3) + :type b: list + + :return: 3a3 matrix of the direction cosine + :rtype: numpy.array + + :Example: + + >>> a = [ + ... [1., 0., 0.], + ... [0., 1., 0.], + ... [0., 0., 1.] + ... ] + >>> b = [ + ... [], + ... [], + ... [] + ... ] + >>> utilities.calcCab(a, b) + ''' + a = np.array(a) + b = np.array(b) + Cab = np.zeros((3, 3)) + for i in range(3): + for j in range(3): + Cab[i][j] = np.dot(a[i], b[j]) + return Cab + + +def distance(p1, p2): + s2 = (p1[0] - p2[0])**2 + (p1[1] - p2[1])**2 + return math.sqrt(s2) + + +def ss(v1, v2, scale=None): + if scale is None: + scale = [1., ] * len(v1) + assert len(scale) == len(v1) + sumsqr = 0.0 + for i in range(len(v1)): + sumsqr = sumsqr + (v1[i] * scale[i] - v2[i] * scale[i])**2 + return sumsqr + + +def parseLayupCode(s_code): + l_code = [] + ilsb = s_code.find('[') + irsb = s_code.find(']') + # print ilsb, irsb + angles = s_code[ilsb + 1:irsb] + # print angles + ns = 0 + if irsb + 1 < len(s_code): + ns = s_code[irsb + 1:].strip('s') + if len(ns) == 0: + ns = '1' + ns = int(ns) + # print ns + + # ilrbs = [] + # irrbs = [] + islash_level1 = [] # index + is_in_brackets = False + for i, c in enumerate(angles): + if c == '(': + is_in_brackets = True + # ilrbs.append(i) + if c == ')': + is_in_brackets = False + # irrbs.append(i) + if (c == '/') and (not is_in_brackets): + islash_level1.append(i) + # print ilrbs + # print irrbs + # print islash_level1 + + islash_level1.append(len(angles)) + list_angles = [] + i = 0 + for j in islash_level1: + list_angles.append(angles[i:j]) + i = j + 1 + + # print list_angles + + for i, a in enumerate(list_angles): + temp = a.split(':') + temp[0] = temp[0].strip('(').strip(')') + # print temp[0] + repeat = 1 + if len(temp) > 1: + repeat = int(temp[-1]) + for k in range(repeat): + if '/' in temp[0]: + temp2 = temp[0].split('/') + for l in temp2: + l_code.append(l) + else: + l_code.append(temp[0]) + + l_code = list(map(float, l_code)) + # print l_code + + for s in range(ns): + temp = l_code[:] + temp.reverse() + # print l_code + # print temp + l_code = l_code + temp + + # for angle in l_code: + # print angle + + return l_code + + +def updateXMLElement(parent, tag, value): + """ Update/Create an element with given tag and value. + + :param parent: Parent element + :type parent: xml.etree.ElementTree.Element + + :param tag: Tag of the element + :type tag: string + + :param value: Value of the element + :type value: string + """ + _se = parent.find(tag) + if _se is None: + _se = et.SubElement(parent, tag) + _se.text = value + + +class Polynomial2DSP(object): + def __init__(self, order, coeffs): + self.order = order + self.coeffs = np.asarray(coeffs) + + def __call__(self, x): + basis = [] + for i in range(self.order+1): + for j in range(i+1): + basis.append(x[0]**(i-j) * x[1]**j) + basis = np.asarray(basis) + + return np.dot(self.coeffs, basis) diff --git a/utils.pyc b/utils.pyc new file mode 100644 index 0000000..77e1dc0 Binary files /dev/null and b/utils.pyc differ