diff --git a/vfo/classTags.py b/vfo/classTags.py new file mode 100644 index 0000000..5d0398c --- /dev/null +++ b/vfo/classTags.py @@ -0,0 +1,10 @@ + + +###### +## Element class tags + +fourNodeEleTags = [31, 32, 52, 53, 55, 156, 157, 203] +triNodeEleTags = [33, 167, 168, 204] +eightNodeBrickEleTags = [36, 38, 43, 44, 56] +MVLEMEleTags = [162, 164, 212, 213] +tetEleTags = [179] diff --git a/vfo/internal_database_functions.py b/vfo/internal_database_functions.py index 3241c58..f6a132d 100644 --- a/vfo/internal_database_functions.py +++ b/vfo/internal_database_functions.py @@ -31,7 +31,9 @@ def _getNodesandElements(): # Check Number of dimensions and intialize variables ndm = len(ops.nodeCoord(nodeList[0])) Nnodes = len(nodeList) + Nele = len(eleList) nodes = np.zeros([Nnodes, ndm + 1]) + eleClassTags = np.zeros([Nele, 2]) # Get Node list for ii, node in enumerate(nodeList): @@ -40,6 +42,7 @@ def _getNodesandElements(): Nele = len(eleList) elements = [None]*Nele + # Generate the element list by looping through all emenemts for ii, ele in enumerate(eleList): @@ -51,9 +54,14 @@ def _getNodesandElements(): tempEle[0] = int(ele) tempEle[1:] = tempNodes - elements[ii] = tempEle + elements[ii] = tempEle + + # Generate element class tags by looping through all elements + for ii, ele in enumerate(eleList): + eleClassTags[ii,0] = ele + eleClassTags[ii,1] = ops.getEleClassTags(ele)[0] - return nodes, elements + return nodes, elements, eleClassTags def _saveNodesandElements(ModelName): """ @@ -81,6 +89,7 @@ def _saveNodesandElements(ModelName): # Consider making these optional arguements nodeName = 'Nodes' eleName = 'Elements' + eleClassName = 'EleClassTags' delim = ' ' fmt = '%.5e' ftype = '.out' @@ -88,7 +97,7 @@ def _saveNodesandElements(ModelName): ODBdir = ModelName+"_ODB" # ODB Dir name # Read noades and elements - nodes, elements = _getNodesandElements() + nodes, elements, eleClassTags = _getNodesandElements() # Sort through the element arrays ele2Node = np.array([ele for ele in elements if len(ele) == 3]) @@ -103,6 +112,8 @@ def _saveNodesandElements(ModelName): ele3File = os.path.join(ODBdir, eleName + "_3Node" + ftype) ele4File = os.path.join(ODBdir, eleName + "_4Node" + ftype) ele8File = os.path.join(ODBdir, eleName + "_8Node" + ftype) + + eleClassTagsFile = os.path.join(ODBdir, eleClassName + ftype) # SaveNodes np.savetxt(nodeFile, nodes, delimiter = delim, fmt = fmt) @@ -112,6 +123,9 @@ def _saveNodesandElements(ModelName): np.savetxt(ele3File, ele3Node, delimiter = delim, fmt = fmt) np.savetxt(ele4File, ele4Node, delimiter = delim, fmt = fmt) np.savetxt(ele8File, ele8Node, delimiter = delim, fmt = fmt) + + # Save Element Class Tags + np.savetxt(eleClassTagsFile, eleClassTags, delimiter = delim, fmt = fmt) def _readNodesandElements(ModelName): @@ -144,6 +158,7 @@ def _readNodesandElements(ModelName): # Consider making these optional arguements nodeName = 'Nodes' eleName = 'Elements' + eleClassName = 'EleClassTags' delim = ' ' dtype ='float32' ftype = '.out' @@ -162,6 +177,8 @@ def _readNodesandElements(ModelName): ele8File = os.path.join(ODBdir, eleName + "_8Node" + ftype) eleFileNames = [ele2File, ele3File, ele4File, ele8File] + + eleClassTagsFile = os.path.join(ODBdir, eleClassName + ftype) ## Load Node information try: @@ -187,11 +204,17 @@ def _readNodesandElements(ModelName): # define the final element array elements = [*TempEle[0],*TempEle[1],*TempEle[2],*TempEle[3]] + ## Load Element Class Tags information + try: + eleClassTags = np.loadtxt(eleClassTagsFile, dtype, delimiter = delim, unpack=False) + except: + print("No element class tag information was found") + # Check if any files were read if elements is []: - raise Exception('No files were found!') + raise Exception('No element information files were found!') - return nodes, elements + return nodes, elements, eleClassTags ################ ModeShapes ############################# diff --git a/vfo/internal_plotting_functions.py b/vfo/internal_plotting_functions.py index 1727e7c..fe96d65 100644 --- a/vfo/internal_plotting_functions.py +++ b/vfo/internal_plotting_functions.py @@ -1,5 +1,6 @@ import numpy as np from mpl_toolkits.mplot3d import Axes3D +from mpl_toolkits.mplot3d.art3d import Poly3DCollection import matplotlib.pyplot as plt import vfo.internal_database_functions as idbf @@ -156,6 +157,88 @@ def _plotQuad3D(iNode, jNode, kNode, lNode, ax, show_element_tags, element, eleS +def _plotTri3D(iNode, jNode, kNode, ax, show_element_tags, element, eleStyle, fillSurface): + ## procedure to render a 2D three node shell element. use eleStyle = "wire" for a wire frame, and "solid" for solid element lines. + ## USe fillSurface = "yes" for color fill in the elements. fillSurface="no" for wireframe. + + ## Initialize varables for matplotlib objects + tempLines = [None] + tempSurface = [None] + tempTag = [None] + + tempLines, = plt.plot((iNode[0], jNode[0], kNode[0], iNode[0]), + (iNode[1], jNode[1], kNode[1], iNode[1]), + (iNode[2], jNode[2], kNode[2], iNode[2]), marker='') + + # update style + if eleStyle == "wire": + plt.setp(tempLines,**WireEle_style) + else: + plt.setp(tempLines,**ele_style) + + # x = [iNode[0],jNode[0],kNode[0]] + # y = [iNode[1],jNode[1],kNode[1]] + # z = [iNode[2],jNode[2],kNode[2]] + # verts = [list(zip(x,y,z))] + # ax.add_collection3d(Poly3DCollection(verts, facecolor='g', alpha=.25)) + + if fillSurface == 'yes': + tempSurface = ax.plot_surface(np.array([[iNode[0], kNode[0]], [jNode[0], kNode[0]]]), + np.array([[iNode[1], kNode[1]], [jNode[1], kNode[1]]]), + np.array([[iNode[2], kNode[2]], [jNode[2], kNode[2]]]), color='g', alpha=.6) + + if show_element_tags == 'yes': + tempTag = ax.text((iNode[0] + jNode[0] + kNode[0])*1.0/3, (iNode[1]+jNode[1]+kNode[1])*1.0/3, + (iNode[2]+jNode[2]+kNode[2])*1.0/3, str(element), **ele_text_style) #label elements + return tempLines, tempSurface, tempTag + + +def _plotTetSurf(nodeCords, ax, fillSurface, eleStyle): + ## This procedure is called by the plotCubeVol() command + aNode = nodeCords[0] + bNode = nodeCords[1] + cNode = nodeCords[2] + + ## Use arrays for less memory and fast code + surfXarray = np.array([[aNode[0], cNode[0]], [bNode[0], cNode[0]]]) + surfYarray = np.array([[aNode[1], cNode[1]], [bNode[1], cNode[1]]]) + surfZarray = np.array([[aNode[2], cNode[2]], [bNode[2], cNode[2]]]) + + ## Initialize varables for matplotlib objects + tempSurface = [None] + + if fillSurface == 'yes': + tempSurface = ax.plot_surface(surfXarray, surfYarray, surfZarray, edgecolor='k', color='g', alpha=.5) + + del aNode, bNode, cNode, surfXarray, surfYarray, surfZarray + + return tempSurface + + +def _plotTetVol(iNode, jNode, kNode, lNode, ax, show_element_tags, element, eleStyle, fillSurface): + ## procedure to render a cubic element, use eleStyle = "wire" for a wire frame, and "solid" for solid element lines. + ## USe fillSurface = "yes" for color fill in the elements. fillSurface="no" for wireframe. + + tempSurfaces = 4*[None] + tempTag = [None] + + # 2D Planer four-node shell elements + # [iNode, jNode, kNode, lNode,iiNode, jjNode, kkNode, llNode] = [*nodesCords] + + tempSurfaces[0] = _plotTetSurf([iNode, jNode, kNode], ax, fillSurface, eleStyle) + tempSurfaces[1] = _plotTetSurf([iNode, jNode, lNode], ax, fillSurface, eleStyle) + tempSurfaces[2] = _plotTetSurf([iNode, kNode, lNode], ax, fillSurface, eleStyle) + tempSurfaces[3] = _plotTetSurf([jNode, kNode, lNode], ax, fillSurface, eleStyle) + + if show_element_tags == 'yes': + tempTag = ax.text((iNode[0]+jNode[0]+kNode[0]+lNode[0])/4, + (iNode[1]+jNode[1]+kNode[1]+lNode[1])/4, + (iNode[2]+jNode[2]+kNode[2]+lNode[2])/4, + str(element), **ele_text_style) #label elements + + return tempSurfaces, tempTag + + def _checkEleLength2D(iNode,jNode): eleLengthCheck = "ELE" diff --git a/vfo/vfo.py b/vfo/vfo.py index a247161..7aba261 100644 --- a/vfo/vfo.py +++ b/vfo/vfo.py @@ -24,6 +24,7 @@ pass from mpl_toolkits.mplot3d import Axes3D +from mpl_toolkits.mplot3d.art3d import Poly3DCollection from math import asin import matplotlib.pyplot as plt import numpy as np @@ -33,6 +34,7 @@ #### CHANGE THESE BEFORE COMMITTING, call them before calling openseespy here #### +import vfo.classTags as classTags import vfo.internal_database_functions as idbf import vfo.internal_plotting_functions as ipltf import openseespy.opensees as ops @@ -278,18 +280,22 @@ def plot_model(*argv,Model="none"): if Model == "none": print("No Model_ODB specified, trying to get data from the active model.") try: - nodeArray, elementArray = idbf._getNodesandElements() + nodeArray, elementArray, eleClassTags = idbf._getNodesandElements() except: raise Exception("No Model_ODB specified. No active model found.") else: print("Reading data from the "+Model+"_ODB.") try: - nodeArray, elementArray = idbf._readNodesandElements(Model) + nodeArray, elementArray, eleClassTags = idbf._readNodesandElements(Model) ## Add eleClassTags to the Tcl data except: raise Exception("No Model_ODB found. No active model found.") nodetags = nodeArray[:,0] + # print(nodeArray) + # print(elementArray) + # print(eleClassTags) + # print(eleClassTags[:,1]) def nodecoords(nodetag): """ @@ -297,6 +303,13 @@ def nodecoords(nodetag): """ i, = np.where(nodeArray[:,0] == float(nodetag)) return nodeArray[int(i),1:] + + def _eleClassTag(eleTag): + """ + Returns element class tag + """ + i, = np.where(eleClassTags[:,0] == float(eleTag)) + return eleClassTags[i,1] # Check if the model is 2D or 3D if len(nodecoords(nodetags[0])) == 2: @@ -324,11 +337,19 @@ def nodecoords(nodetag): ipltf._plotTri2D(iNode, jNode, kNode, ax, show_element_tags, eleTag, ele_style, fillSurface='yes') if len(Nodes) == 4: - # 2D Planer four-node shell elements - iNode = nodecoords(Nodes[0]) - jNode = nodecoords(Nodes[1]) - kNode = nodecoords(Nodes[2]) - lNode = nodecoords(Nodes[3]) + if _eleClassTag(eleTag) in classTags.MVLEMEleTags: + # 2D Planer four-node shell elements + iNode = nodecoords(Nodes[0]) + jNode = nodecoords(Nodes[1]) + kNode = nodecoords(Nodes[3]) + lNode = nodecoords(Nodes[2]) + + else: + ## 2D Planer four-node shell elements + iNode = nodecoords(Nodes[0]) + jNode = nodecoords(Nodes[1]) + kNode = nodecoords(Nodes[2]) + lNode = nodecoords(Nodes[3]) ipltf._plotQuad2D(iNode, jNode, kNode, lNode, ax, show_element_tags, eleTag, ele_style, fillSurface='yes') @@ -359,14 +380,43 @@ def nodecoords(nodetag): ipltf._plotBeam3D(iNode, jNode, ax, show_element_tags, eleTag, "solid") - if len(Nodes) == 4: - # 3D four-node Quad/shell element + if len(Nodes) == 3: + # 3D Planer three-node shell elements iNode = nodecoords(Nodes[0]) jNode = nodecoords(Nodes[1]) kNode = nodecoords(Nodes[2]) - lNode = nodecoords(Nodes[3]) - ipltf._plotQuad3D(iNode, jNode, kNode, lNode, ax, show_element_tags, eleTag, ele_style, fillSurface='yes') + ipltf._plotTri3D(iNode, jNode, kNode, ax, show_element_tags, eleTag, ele_style, fillSurface='yes') + + if len(Nodes) == 4: + if _eleClassTag(eleTag) in classTags.MVLEMEleTags: + ## 3D four-node Quad/shell element + ## Reverse the node connectivity for MVLEM3D elements + iNode = nodecoords(Nodes[0]) + jNode = nodecoords(Nodes[1]) + kNode = nodecoords(Nodes[3]) + lNode = nodecoords(Nodes[2]) + + ipltf._plotQuad3D(iNode, jNode, kNode, lNode, ax, show_element_tags, eleTag, ele_style, fillSurface='yes') + + elif _eleClassTag(eleTag) in classTags.fourNodeEleTags: + ## 3D Planer four-node shell elements + iNode = nodecoords(Nodes[0]) + jNode = nodecoords(Nodes[1]) + kNode = nodecoords(Nodes[2]) + lNode = nodecoords(Nodes[3]) + + ipltf._plotQuad3D(iNode, jNode, kNode, lNode, ax, show_element_tags, eleTag, ele_style, fillSurface='yes') + + elif _eleClassTag(eleTag) in classTags.tetEleTags: + ## If the element is a four-node tetrahedron + iNode = nodecoords(Nodes[0]) + jNode = nodecoords(Nodes[1]) + kNode = nodecoords(Nodes[2]) + lNode = nodecoords(Nodes[3]) + + ipltf._plotTetVol(iNode, jNode, kNode, lNode, ax, show_element_tags, eleTag, 'solid', fillSurface='yes') + if len(Nodes) == 8: # 3D eight-node Brick element diff --git a/vfo_test.py b/vfo_test.py new file mode 100644 index 0000000..845f30f --- /dev/null +++ b/vfo_test.py @@ -0,0 +1,322 @@ + + +#################################################### +## +## Test Examples +## +#################################################### + +import numpy as np +import openseespy.opensees as ops +import vfo.vfo as vfo + +import math + +def portalframe2d(): + + print("portal frame 2D test") + + ops.wipe() + + # set modelbuilder + ops.model('basic', '-ndm', 2, '-ndf', 3) + + ############################################ + ### Units and Constants ################### + ############################################ + + inch = 1; + kip = 1; + sec = 1; + + # Dependent units + sq_in = inch*inch; + ksi = kip/sq_in; + ft = 12*inch; + + # Constants + g = 386.2*inch/(sec*sec); + pi = math.acos(-1); + + ####################################### + ##### Dimensions + ####################################### + + # Dimensions Input + H_story=10.0*ft; + W_bayX=16.0*ft; + W_bayY_ab=5.0*ft+10.0*inch; + W_bayY_bc=8.0*ft+4.0*inch; + W_bayY_cd=5.0*ft+10.0*inch; + + # Calculated dimensions + W_structure=W_bayY_ab+W_bayY_bc+W_bayY_cd; + + # ########################### + # ##### Nodes + # ########################### + + # Create All main nodes + ops.node(1, 0.0, 0.0) + ops.node(2, W_bayX, 0.0) + + ops.node(11, 0.0, H_story) + ops.node(12, W_bayX, H_story) + + # ############### + # Constraints + # ############### + + ops.fix(1, 1, 1, 1) + ops.fix(2, 1, 1, 1) + + # ####################### + # ### Elements + # ####################### + + ColTransfTag=1 + BeamTranfTag=2 + + ops.geomTransf('Linear', ColTransfTag) + ops.geomTransf('Linear', BeamTranfTag) + + # Assign Elements ############## + # ## Add non-linear column elements + ops.element('elasticBeamColumn', 1, 1, 11, 20., 1000., 1000., ColTransfTag, '-mass', 0.0) + ops.element('elasticBeamColumn', 2, 2, 12, 20., 1000., 1000., ColTransfTag, '-mass', 0.0) + + # ### Add linear main beam elements, along x-axis + ops.element('elasticBeamColumn', 101, 11, 12, 20., 1000., 1000., BeamTranfTag, '-mass', 0.0) + + # Visualize the model + return vfo.plot_model() + + +def portalframe3d(): + print("portal frame 3D test") + + ops.wipe() + + # set modelbuilder + ops.model('basic', '-ndm', 3, '-ndf', 6) + + ############################################ + ### Units and Constants ################### + ############################################ + + inch = 1; + kip = 1; + sec = 1; + + # Dependent units + sq_in = inch*inch; + ksi = kip/sq_in; + ft = 12*inch; + + # Constants + g = 386.2*inch/(sec*sec); + pi = math.acos(-1); + + ####################################### + ##### Dimensions + ####################################### + + # Dimensions Input + H_story=10.0*ft; + W_bayX=16.0*ft; + W_bayY_ab=5.0*ft+10.0*inch; + W_bayY_bc=8.0*ft+4.0*inch; + W_bayY_cd=5.0*ft+10.0*inch; + + # Calculated dimensions + W_structure=W_bayY_ab+W_bayY_bc+W_bayY_cd; + + # ########################### + # ##### Nodes + # ########################### + + # Create All main nodes + ops.node(1, 0.0, 0.0, 0.0) + ops.node(2, W_bayX, 0.0, 0.0) + + ops.node(11, 0.0, H_story, 0.0) + ops.node(12, W_bayX, H_story, 0.0) + + # ############### + # Constraints + # ############### + + ops.fix(1, 1, 1, 1, 1, 1, 1) + ops.fix(2, 1, 1, 1, 1, 1, 1) + + # ####################### + # ### Elements + # ####################### + + ColTransfTag=1 + BeamTranfTag=2 + + ops.geomTransf('Linear', ColTransfTag, 1, 0, 0) + ops.geomTransf('Linear', BeamTranfTag, 0, 0, 1) + + # Assign Elements ############## + # ## Add non-linear column elements + ops.element('elasticBeamColumn', 1, 1, 11, 20., 1000., 1000., 1000., 1000., 1000., ColTransfTag, '-mass', 0.0) + ops.element('elasticBeamColumn', 2, 2, 12, 20., 1000., 1000., 1000., 1000., 1000.,ColTransfTag, '-mass', 0.0) + + # ### Add linear main beam elements, along x-axis + ops.element('elasticBeamColumn', 101, 11, 12, 20., 1000., 1000., 1000., 1000., 1000.,BeamTranfTag, '-mass', 0.0) + + # Visualize the model + return vfo.plot_model() + + +def tri3d(): + print("portal frame 3D test") + + ops.wipe() + + # set modelbuilder + ops.model('basic', '-ndm', 3, '-ndf', 3) + + # Create All main nodes + ops.node(1, 0.0, 0.0, 0.0) + ops.node(2, 10.0, 0.0, 0.0) + ops.node(3, 0.0, 10.0, 0.0) + + # ############### + # Constraints + # ############### + + ops.fix(1, 1, 1, 1) + ops.fix(2, 1, 1, 1) + + ops.nDMaterial("ElasticIsotropic", 1, 1000.0, 0.25, 6.75) + + # Assign Elements ############## + # ## Add non-linear column elements + ops.element('Tri31', 1, 1,2,3, 0.1, 'PlaneStress', 1) + + # Visualize the model + return vfo.plot_model() + + +def quad2d(): + print("quad3d 3D test") + + ops.wipe() + + # set modelbuilder + ops.model('basic', '-ndm', 2, '-ndf', 2) + + # ########################### + # ##### Nodes + # ########################### + + # Create All main nodes + ops.node(1, 0.0, 0.0) + ops.node(2, 50.0, 0.0) + ops.node(3, 0.0, 30.0) + ops.node(4, 50.0, 30.0) + + # ops.fix(1, 1, 1, 1, 1, 1, 1) + # ops.fix(2, 1, 1, 1, 1, 1, 1) + + ops.fix(1, 1, 1) + ops.fix(2, 1, 1) + + ops.nDMaterial("ElasticIsotropic", 1, 1000.0, 0.25, 6.75) + + # Assign Elements ############## + # ## Add non-linear column elements + ops.element('quad', 1, 1,2,4,3, 0.1, 'PlaneStress', 1) + + # vfo.createODB("Test_mvlem3d", "none", 0) + # Visualize the model + return vfo.plot_model() + + + +def mvlem_3d(): + print("MVLEM 3D test") + + ops.wipe() + + # set modelbuilder + ops.model('basic', '-ndm', 3, '-ndf', 6) + + # ########################### + # ##### Nodes + # ########################### + + # Create All main nodes + ops.node(1, 0.0, 0.0, 0.0) + ops.node(2, 50.0, 0.0, 0.0) + ops.node(3, 0.0, 0.0, 30.0) + ops.node(4, 50.0, 0.0, 30.0) + ops.node(5, 0.0, 0.0, 60.0) + ops.node(6, 50.0, 0.0, 60.0) + + ops.fix(1, 1, 1, 1, 1, 1, 1) + ops.fix(2, 1, 1, 1, 1, 1, 1) + + ops.uniaxialMaterial('Elastic', 1, 314705) + # Concrete Materials + ops.uniaxialMaterial('Concrete02', 201, -7.934, -0.0023, 0, -0.01, 0.079, 0.356292015066294, 253.858060734734) + + # Steel Materials + ops.uniaxialMaterial('SteelMPF', 301, 68.313, 68.313, 27847.246, 0.0055, 0.0055, 20, 0.925, 0.15) + + # Assign Elements ############## + ops.element('MVLEM_3D', 1,1,2,4,3, 10, '-thick', *[3.937, 3.937, 3.937, 3.937, 3.937, 3.937, 3.937, 3.937, 3.937, 3.937], + '-width', *[5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0], + '-rho', *[0.0387, 0.0387, 0.00342, 0.00342, 0.00342, 0.00342, 0.00342, 0.00342, 0.0226, 0.0226], + '-matConcrete', *[201, 201, 201, 201, 201, 201, 201, 201, 201, 201], + '-matSteel', *[301, 301, 301, 301, 301, 301, 301, 301, 301, 301], '-matShear', 1) + ops.element('MVLEM_3D', 2,3,4,6,5, 10, '-thick', *[3.937, 3.937, 3.937, 3.937, 3.937, 3.937, 3.937, 3.937, 3.937, 3.937], + '-width', *[5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0], + '-rho', *[0.0387, 0.0387, 0.00342, 0.00342, 0.00342, 0.00342, 0.00342, 0.00342, 0.0226, 0.0226], + '-matConcrete', *[201, 201, 201, 201, 201, 201, 201, 201, 201, 201], + '-matSteel', *[301, 301, 301, 301, 301, 301, 301, 301, 301, 301], '-matShear', 1) + + # vfo.createODB("Test_mvlem3d", "none", 0) + # Visualize the model + return vfo.plot_model() + +def tetra(): + print("Tetrahedral test") + + ops.wipe() + + # set modelbuilder + ops.model('basic', '-ndm', 3, '-ndf', 3) + + # Create All main nodes + ops.node(1, 0.0, 0.0, 0.0) + ops.node(2, 10.0, 0.0, 0.0) + ops.node(3, 5.0, 10.0, 0.0) + ops.node(4, 5.0, 5.0, 10.0) + + # ############### + # Constraints + # ############### + + ops.fix(1, 1, 1, 1) + ops.fix(2, 1, 1, 1) + + ops.nDMaterial("ElasticIsotropic", 1, 1000.0, 0.25, 6.75) + + # Assign Elements ############## + # ## Add non-linear column elements + ops.element('FourNodeTetrahedron', 1, 1,2,3,4, 1) + + # Visualize the model + return vfo.plot_model() + + +fig1 = portalframe2d() +fig2 = portalframe3d() +fig3 = quad2d() +fig3 = mvlem_3d() +fog4 = tri3d() +fog5 = tetra()