From b35a19815e50e77540707e601c196e696bad94c3 Mon Sep 17 00:00:00 2001 From: Su Tian Date: Thu, 29 Jul 2021 17:55:40 -0400 Subject: [PATCH] 0729 --- msgpi/sg.py | 540 +++++++++++++++++++++++++++++----------------------- 1 file changed, 305 insertions(+), 235 deletions(-) diff --git a/msgpi/sg.py b/msgpi/sg.py index 049b554..1be789f 100644 --- a/msgpi/sg.py +++ b/msgpi/sg.py @@ -104,241 +104,6 @@ def summary(self): print(row) - def getBeamProperty(self, label): - """Get beam properties using specific labels. - - Parameters - ---------- - label : str - Label of the property that will be returned. - - Mass - - * msijo - Entry (i, j) of the 6x6 mass matrix at the origin - * msijc - Entry (i, j) of the 6x6 mass matrix at the mass center - * mpl - Mass per unit length - * mmoi1/mmoi2/mmoi3 - Mass moment of inertia about x1/x2/x3 axis - - Stiffness - - * stfijc - Entry (i, j) of the 4x4 classical stiffness matrix - * stfijr - Entry (i, j) of the 6x6 refined stiffness matrix - * eac/ear - Axial stiffness of the classical/refined model - * gjc/gjr - Torsional stiffness of the classical/refined model - * ei2c/eifc/ei2r/eifr - Bending stiffness around x2 (flapwise) of the classical/refined model - * ei3c/eicc/ei3r/eicr - Bending stiffness around x3 (chordwise or lead-lag) of the classical/refined model - - Compliance - - * cmpijc - Entry (i, j) of the 4x4 classical compliance matrix - * cmpijr - Entry (i, j) of the 6x6 refined compliance matrix - - Centers - - * mcy/mc2 - y (or x2) component of the mass center - * mcz/mc3 - z (or x3) component of the mass center - * tcy/tc2 - y (or x2) component of the tension center - * tcz/tc3 - z (or x3) component of the tension center - * scy/sc2 - y (or x2) component of the shear center - * scz/sc3 - z (or x3) component of the shear center - """ - - mm = self.mass_origin - stf_c = self.stiffness - cmp_c = self.compliance - stf_r = self.stiffness_refined - cmp_r = self.compliance_refined - - # if type(center).__name__ == 'str': - # if center == 'tc': - # mm, stf_c, cmp_c, stf_r, cmp_r = calcOffsetBeamProperty(self.tension_center[1], self.tension_center[2]) - # elif center == 'sc': - # mm, stf_c, cmp_c, stf_r, cmp_r = calcOffsetBeamProperty(self.shear_center[1], self.shear_center[2]) - # elif center == 'mc': - # mm, stf_c, cmp_c, stf_r, cmp_r = calcOffsetBeamProperty(self.mass_center[1], self.mass_center[2]) - - # Mass - if label.startswith('ms'): - return mm[int(label[2])-1][int(label[3])-1] - if label == 'mpl': - return self.density - if label == 'mmoi1': - return mm[3, 3] - if label == 'mmoi2': - return mm[4, 4] - if label == 'mmoi3': - return mm[5, 5] - - # Stiffness - if label.startswith('stf'): - if label[-1] == 'c': - return stf_c[int(label[3])-1][int(label[4])-1] - elif label[-1] == 'r': - return stf_r[int(label[3])-1][int(label[4])-1] - if label.startswith('ea'): - # Axial stiffness - if label[-1] == 'c': - return stf_c[0][0] - elif label[-1] == 'r': - return stf_r[0][0] - if label.startswith('gj'): - # Torsional stiffness - if label[-1] == 'c': - return stf_c[1][1] - elif label[-1] == 'r': - return stf_r[3][3] - if label.startswith('ei2') or label.startswith('eif'): - # Bending stiffness around x2 (flapwise) - if label[-1] == 'c': - return stf_c[2][2] - elif label[-1] == 'r': - return stf_r[4][4] - if label.startswith('ei3') or label.startswith('eic'): - # Bending stiffness around x3 (chordwise or lead-lag) - if label[-1] == 'c': - return stf_c[3][3] - elif label[-1] == 'r': - return stf_r[5][5] - - # Compliance - if label.startswith('cmp'): - if label[-1] == 'c': - return cmp_c[int(label[3])-1][int(label[4])-1] - elif label[-1] == 'r': - return cmp_r[int(label[3])-1][int(label[4])-1] - - # Various centers - if label == 'mcy' or label == 'mc2': - return self.mass_center[1] - if label == 'mcz' or label == 'mc3': - return self.mass_center[2] - if label == 'tcy' or label == 'tc2': - return self.tension_center[1] - if label == 'tcz' or label == 'tc3': - return self.tension_center[2] - if label == 'scy' or label == 'sc2': - return self.shear_center[1] - if label == 'scz' or label == 'sc3': - return self.shear_center[2] - - - def calcOffsetBeamProperty(self, offset_x2=0.0, offset_x3=0.0): - - # Offset mass matrix - mm_o = np.asarray(self.mass_mc) - if (offset_x2 != self.mass_center[1]) or (offset_x3 != self.mass_center[2]): - # mm_c = np.asarray(self.mass_mc) - mu = mm_o[0, 0] - mi_c = mm_o[3:, 3:] - - x2 = self.mass_center[1] - offset_x2 - x3 = self.mass_center[2] - offset_x3 - r_tilde = np.array([ - [0, -x3, x2], - [x3, 0, 0], - [-x2, 0, 0] - ]) - - mm_o[:3, 3:] = mu * r_tilde.T - mm_o[3:, :3] = mu * r_tilde - - # I_o = I_c + m * r_tilde.r_tilde^T - mm_o[3:, 3:] = mm_o[3:, 3:] + mu * np.dot(r_tilde, r_tilde.T) - - - # Offset stiffness and compliance - trfm_4 = np.eye(4) - trfm_6 = np.eye(6) - - trfm_4[2, 0] = offset_x3 - trfm_4[3, 0] = -offset_x2 - - trfm_6[4, 0] = offset_x3 - trfm_6[5, 0] = -offset_x2 - trfm_6[3, 1] = -offset_x3 - trfm_6[3, 2] = offset_x2 - - cmp_4 = np.asarray(self.compliance) - cmp_6 = np.asarray(self.compliance_refined) - - cmp_4_trf = np.dot(trfm_4.T, np.dot(cmp_4, trfm_4)) - cmp_6_trf = np.dot(trfm_6.T, np.dot(cmp_6, trfm_6)) - - stf_4_trf = np.linalg.inv(cmp_4_trf) - stf_6_trf = np.linalg.inv(cmp_6_trf) - - return mm_o, stf_4_trf, cmp_4_trf, stf_6_trf, cmp_6_trf - - - def offsetBeamRefCenter(self, offset_x2, offset_x3): - """Offset the beam reference center and recalculate beam properties. - - Parameters - ---------- - offset_x2 : float - x2 of the offset of the new center with respect to the current one. - offset_x3 : float - x3 of the offset of the new center with respect to the current one. - - """ - - # offset = [0.0, offset[0], offset[1]] - - # Offset mass matrix - mm_o = np.asarray(self.mass_mc) - if (offset_x2 != self.mass_center[1]) or (offset_x3 != self.mass_center[2]): - # mm_c = np.asarray(self.mass_mc) - mu = mm_o[0, 0] - mi_c = mm_o[3:, 3:] - - x2 = self.mass_center[1] - offset_x2 - x3 = self.mass_center[2] - offset_x3 - r_tilde = np.array([ - [0, -x3, x2], - [x3, 0, 0], - [-x2, 0, 0] - ]) - - mm_o[:3, 3:] = mu * r_tilde.T - mm_o[3:, :3] = mu * r_tilde - - # I_o = I_c + m * r_tilde.r_tilde^T - mm_o[3:, 3:] = mm_o[3:, 3:] + mu * np.dot(r_tilde, r_tilde.T) - self.mass_origin = mm_o - self.mass_center[1] -= offset_x2 - self.mass_center[2] -= offset_x3 - - - # Offset stiffness and compliance - trfm_4 = np.eye(4) - trfm_6 = np.eye(6) - - trfm_4[2, 0] = offset_x3 - trfm_4[3, 0] = -offset_x2 - - trfm_6[4, 0] = offset_x3 - trfm_6[5, 0] = -offset_x2 - trfm_6[3, 1] = -offset_x3 - trfm_6[3, 2] = offset_x2 - - cmp_4 = np.asarray(self.compliance) - cmp_6 = np.asarray(self.compliance_refined) - - self.compliance = np.dot(trfm_4.T, np.dot(cmp_4, trfm_4)) - self.compliance_refined = np.dot(trfm_6.T, np.dot(cmp_6, trfm_6)) - - self.stiffness = np.linalg.inv(self.compliance) - self.stiffness_refined = np.linalg.inv(self.compliance_refined) - - self.tension_center[1] -= offset_x2 - self.tension_center[2] -= offset_x3 - - self.shear_center[1] -= offset_x2 - self.shear_center[2] -= offset_x3 - - return - - @@ -666,6 +431,48 @@ def offsetBeamRefCenter(self, offset_x2, offset_x3): +class ShellProperty(MaterialSection): + """ + """ + + + + + + + + + +class MaterialProperty(MaterialSection): + """ + """ + def __init__(self, name=''): + MaterialSection.__init__(self, 3) + self.name = name + self.e1 = None + self.e2 = None + self.e3 = None + self.g12 = None + self.g13 = None + self.g23 = None + self.nu12 = None + self.nu13 = None + self.nu23 = None + + def assignConstants(self, consts): + if len(consts) == 9: + self.e1, self.e2, self.e3 = consts[:3] + self.g12, self.g13, self.g23 = consts[3:6] + self.nu12, self.nu13, self.nu23 = consts[6:] + + + + + + + + + class StructureGene(object): """ A finite element level structure gene model in the theory of MSG. @@ -923,3 +730,266 @@ def writeGmshMsh(self, fo, nid_begin=1, eid_begin=1, loc=[0, 0, 0], *args, **kwa return + + def __writeSCNodes(self, fobj, sfi, sff): + if self.sgdim == 1: + nid = self.elements[self.elementids1d[0]][0] + fobj.write(sfi.format(nid)) + ioutil.writeFormatFloats(fobj, self.nodes[nid]) + for eid in self.elementids1d: + if len(self.elements[eid]) > 2: + # node 3 + nid = self.elements[eid][2] + fobj.write(sfi.format(nid)) + ioutil.writeFormatFloats(fobj, self.nodes[nid]) + if len(self.elements[eid]) == 5: + # node 5 + nid = self.elements[eid][4] + fobj.write(sfi.format(nid)) + ioutil.writeFormatFloats(fobj, self.nodes[nid]) + if len(self.elements[eid]) > 3: + # node 4 + nid = self.elements[eid][3] + fobj.write(sfi.format(nid)) + ioutil.writeFormatFloats(fobj, self.nodes[nid]) + # node 2 + nid = self.elements[eid][1] + fobj.write(sfi.format(nid)) + ioutil.writeFormatFloats(fobj, self.nodes[nid]) + else: + for nid, ncoord in self.nodes.items(): + fobj.write(sfi.format(nid)) + ioutil.writeFormatFloats(fobj, ncoord) + return + + + + + def __writeSCElements(self, fobj, sfi): + for eid in self.elementids1d: + elem = np.zeros(7, dtype=int) + elem[0] = eid + elem[1] = self.elem_prop[eid] # mid/cid + for i, nid in enumerate(self.elements[eid]): + elem[i+2] = nid + ioutil.writeFormatIntegers(fobj, elem) + + for eid in self.elementids2d: + elem = np.zeros(11, dtype=int) + elem[0] = eid + elem[1] = self.elem_prop[eid] + elem_type = 'quad' + if (len(self.elements[eid]) == 3) or (len(self.elements[eid]) == 6): + elem_type = 'tri' + for i, nid in enumerate(self.elements[eid]): + if (i >= 3) and (elem_type == 'tri'): + elem[i+3] = nid + else: + elem[i+2] = nid + ioutil.writeFormatIntegers(fobj, elem) + + for eid in self.elementids3d: + elem = np.zeros(22, dtype=int) + elem[0] = eid + elem[1] = self.elem_prop[eid] + elem_type = 'hexa' + if (len(self.elements[eid]) == 4) or (len(self.elements[eid]) == 10): + elem_type = 'tetra' + for i, nid in enumerate(self.elements[eid]): + if (i >= 4) and (elem_type == 'tetra'): + elem[i+3] = nid + else: + elem[i+2] = nid + ioutil.writeFormatIntegers(fobj, elem) + + return + + + + + def __writeSCElementOrientations(self, fobj, sfi, sff): + for eid, orient in self.elem_orient.items(): + fobj.write(sfi.format(eid)) + ioutil.writeFormatFloats(fobj, np.hstack(orient)) + return + + + + + def __writeSCMOCombos(self, fobj, sfi, sff): + for cid, combo in self.mocombos.items(): + fobj.write((sfi + sfi + sff + '\n').format(cid, combo[0], combo[1])) + return + + + + + def __writeSCMaterials(self, fobj, sfi, sff): + for mid, m in self.materials.items(): + # itype = 0 + # if m.type == 'orthotropic': + # itype = 1 + # elif m.type == 'anisotropic': + # itype = 2 + mp = m.eff_props[3] + ioutil.writeFormatIntegers(fobj, (mid, mp['type'], 1)) + ioutil.writeFormatFloats(fobj, (0, mp['density'])) + if mp['type'] == 0: + mpc = mp['constants'] + ioutil.writeFormatFloats(fobj, [mpc['e'], mpc['nu']]) + elif mp['type'] == 1: + mpc = mp['constants'] + ioutil.writeFormatFloats(fobj, [mpc['e1'], mpc['e2'], mpc['e3']]) + ioutil.writeFormatFloats(fobj, [mpc['g12'], mpc['g13'], mpc['g23']]) + ioutil.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(self, fn, sfi, sff): + with open(fn, 'w') as fobj: + # Extra inputs for dimensionally reducible structures + if (self.smdim == 1) or (self.smdim == 2): + # model (0: classical, 1: shear refined) + fobj.write((sfi + '\n').format(self.model)) + if self.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 self.smdim == 2: # shell + # initial twist/curvatures + fobj.write((sff * 2 + '\n').format( + self.initial_curvature[0], self.initial_curvature[1] + )) + fobj.write('\n') + + # Head + fobj.write((sfi * 4 + '\n').format( + self.physics, + self.degen_element, + self.trans_element, + self.nonuniform_temperature + )) + fobj.write((sfi * 6 + '\n').format( + self.sgdim, + len(self.nodes), + len(self.elements), + len(self.materials), + self.num_slavenodes, + len(self.mocombos) + )) + fobj.write('\n') + + # Nodal coordinates + self.__writeSCNodes(self, fobj, sfi, sff) + fobj.write('\n') + + # Element connectivities + self.__writeSCElements(self, fobj, sfi) + fobj.write('\n') + + if self.trans_element != 0: + self.__writeSCElementOrientations(self, fobj, sfi, sff) + fobj.write('\n') + + # Material-orientation combinations + if len(self.mocombos) > 0: + self.__writeSCMOCombos(self, fobj, sfi, sff) + fobj.write('\n') + + # Material + self.__writeSCMaterials(self, fobj, sfi, sff) + fobj.write('\n') + + # Omega + fobj.write((sff + '\n').format(self.omega)) + fobj.write('\n') + + fobj.write('\n') + + + + + + + + + + def writeSCInD(self, fn, sfi, sff): + with open(fn, 'w') as fobj: + ioutil.writeFormatFloats(fobj, self.global_displacements, sff[2:-1]) + ioutil.writeFormatFloatsMatrix(fobj, self.global_rotations, sff[2:-1]) + fobj.write((sfi+'\n').format(self.global_loads_type)) + ioutil.writeFormatFloats(fobj, self.global_loads, sff[2:-1]) + + + + + + + + + + def writeSCInF(self, fn, sfi, sff): + with open(fn, 'w') as fobj: + # Addtional material properties + for i, m in self.materials.items(): + ioutil.writeFormatIntegers( + fobj, + (m.strength['criterion'], len(m.strength['constants'])), + sfi[2:-1] + ) + fobj.write((sff+'\n').format(m.strength['chara_len'])) + ioutil.writeFormatFloats(fobj, m.strength['constants'], sff[2:-1]) + + # Load type + fobj.write((sfi+'\n').format(self.global_loads_type)) + + # (fe) Axes + + # (fi) Loads + ioutil.writeFormatFloats(fobj, self.global_loads, sff[2:-1]) + + + + + + + + + + def writeSCIn(self, fn, analysis='h'): + """ Write SG input files for SwiftComp + + :param analysis: Type of analysis. + - 'h': Homogenization + - 'd' or 'l': Dehomogenization + - 'f': Failure analysis + """ + + + # string format + sfi = '{:8d}' + sff = '{:16.6E}' + + if 'h' in analysis: + self.writeSCInH(self, fn, sfi, sff) + elif ('d' in analysis) or ('l' in analysis): + self.writeSCInD(self, fn, sfi, sff) + + return fn + +