Skip to content

Commit

Permalink
Added plot functions
Browse files Browse the repository at this point in the history
-MVLEM
- Tetrahedron
- Tri3D
  • Loading branch information
u-anurag committed Jan 3, 2022
1 parent df2146c commit ae2023d
Show file tree
Hide file tree
Showing 5 changed files with 504 additions and 16 deletions.
10 changes: 10 additions & 0 deletions vfo/classTags.py
Original file line number Diff line number Diff line change
@@ -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]
33 changes: 28 additions & 5 deletions vfo/internal_database_functions.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand All @@ -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):
Expand All @@ -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):
"""
Expand Down Expand Up @@ -81,14 +89,15 @@ def _saveNodesandElements(ModelName):
# Consider making these optional arguements
nodeName = 'Nodes'
eleName = 'Elements'
eleClassName = 'EleClassTags'
delim = ' '
fmt = '%.5e'
ftype = '.out'

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])
Expand All @@ -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)
Expand All @@ -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):
Expand Down Expand Up @@ -144,6 +158,7 @@ def _readNodesandElements(ModelName):
# Consider making these optional arguements
nodeName = 'Nodes'
eleName = 'Elements'
eleClassName = 'EleClassTags'
delim = ' '
dtype ='float32'
ftype = '.out'
Expand All @@ -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:
Expand All @@ -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 #############################
Expand Down
83 changes: 83 additions & 0 deletions vfo/internal_plotting_functions.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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"
Expand Down
72 changes: 61 additions & 11 deletions vfo/vfo.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -278,25 +280,36 @@ 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):
"""
Returns an array of node coordinates: works like nodeCoord() in opensees.
"""
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:
Expand Down Expand Up @@ -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')

Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit ae2023d

Please sign in to comment.