# # This script should only be used for DEM cluster visualization # The script takes an xxxx.odb and generates and new file named xxxx_clst.odb. # The new file xxxx_clst.odb can be opened in Abaqus/Viewer to visualize # the dem clusters in the model. # import sys import odbAccess from abaqusConstants import * from odbAccess import * from odbSection import * from symbolicConstants import * import math import time import numpy as np import itertools linearg = sys.argv[1:] jobName, extension = linearg[0].split('.') odbName = jobName + '.odb' try: session.odbs[odbName].close() except: pass from shutil import copyfile src = odbName dst = jobName + "_clst" + '.odb' startTime = time.time() print("****************\nCopying the Odb... at{}".format(startTime)) try: copyfile(src, dst) except: pass #------------------------------------------------------------------------------- # generate child nodes and elements childParticleRadiusSection = [] childParticleArmSection = [] def createChildParticles(odb, secIndex, section, Rp, childParticles): a=odb.rootAssembly OriName = a.datumCsyses.keys() elset_Name = section.region.name OriNameCls = 'ORI-'+elset_Name lOri = 0 for i in range(len(OriName)): if (OriNameCls == OriName[i]): lOri = 1 vec = [[0 for i in range(3)] for j in range(3)] vec[0][0] = 1.0 vec[0][1] = 0.0 vec[0][2] = 0.0 vec[1][0] = 0.0 vec[1][1] = 1.0 vec[1][2] = 0.0 vec[2][0] = 0.0 vec[2][1] = 0.0 vec[2][2] = 1.0 if lOri == 1: xAxis = a.datumCsyses[OriNameCls].xAxis yAxis = a.datumCsyses[OriNameCls].yAxis zAxis = a.datumCsyses[OriNameCls].zAxis ClstOri(xAxis, yAxis, zAxis,vec) instanceName = section.region.elements[0].instanceName inst = a.instances[instanceName] partName = section.sectionName + '_child' part = odb.Part(partName, embeddedSpace=THREE_D, type=DEFORMABLE_BODY) childParticlesLocation = [] childParticlesRadius = [] childArm = [] for child in childParticles: alpha = float(child[0]) beta = float(child[1]) theta = float(child[2]) phi = float(child[3]) r = alpha * Rp R = beta * Rp x = r * math.sin((phi/180)*math.pi)* math.cos((theta/180)*math.pi) y = r * math.sin((phi/180)*math.pi)* math.sin((theta/180)*math.pi) z = r * math.cos((phi/180)*math.pi) VecX = vec[0][0]*x + vec[0][1]*y + vec[0][2]*z VecY = vec[1][0]*x + vec[1][1]*y + vec[1][2]*z VecZ = vec[2][0]*x + vec[2][1]*y + vec[2][2]*z childParticlesLocation.append((VecX,VecY,VecZ)) childParticlesRadius.append(R) childParticleRadiusSection.append(childParticlesRadius) childParticleArmSection.append(childParticlesLocation) nodeStartLabel = 1; elemStartLabel = 1 nodeLabel = nodeStartLabel elemLabel = elemStartLabel nodeData = []; elementData = [] for elem in section.region.elements: nodeL = elem.connectivity[0] x0,y0,z0 = inst.getNodeFromLabel(nodeL).coordinates for x1,y1,z1 in childParticlesLocation: x = x0+x1 y = y0+y1 z = z0+z1 nodeData.append((nodeLabel, x, y, z)) elementData.append((elemLabel, nodeLabel)) nodeLabel += 1 elemLabel += 1 print('number of children particles created for section %s: %d'%(section.sectionName, len(nodeData))) print('adding children nodes...at', time.time()) part.addNodes(nodeData=nodeData,nodeSetName= 'nset') print('adding children elems...at', time.time()) part.addElements(elementData=elementData, type='PD3D', elementSetName='elset') newInst=a.Instance(name= partName + '-1', object=part) odb.save() #------------------------------------------------------------------------------- # compute the current configuration of clusters def applyDisplacement(odb, childInstanceName, parentSection, ChildParticleRad, childParticleArm): a = odb.rootAssembly nodeStartLabel = 1 inst = a.instances[childInstanceName] numChildElems = len(inst.nodes) parentRegion = parentSection.region parentElems = parentRegion.elements numParentElems = len(parentElems) numChildPerParent = int(numChildElems/numParentElems) nodeLabel = nodeStartLabel iscField = False parentInstanceName = parentElems[0].instanceName parentNodes = a.instances[parentInstanceName].nodeSets[parentRegion.name] for step in odb.steps.values(): for frameNumber in range(len(step.frames)): #len(step.frames) dispData = [] pradData = [] coordData = [] nodeLabels = [] uField = step.frames[frameNumber].fieldOutputs['U'] uFieldSub = uField.getSubset(region = parentNodes) uRField = step.frames[frameNumber].fieldOutputs['UR'] uRFieldSub = uRField.getSubset(region = parentNodes) iscField = False if('COORD' in step.frames[frameNumber].fieldOutputs): iscField = True if(iscField): cField = step.frames[frameNumber].fieldOutputs['COORD'] cFieldSub = cField.getSubset(region = parentNodes) pradField = step.frames[frameNumber].fieldOutputs['PRAD'] pradFieldSub = pradField.getSubset(region = parentNodes) if frameNumber == 0: #populate the nodeLabels just for the first frame nodeLabel = nodeStartLabel if(iscField): for item in zip(uFieldSub.values, uRFieldSub.values, cFieldSub.values): disp = tuple(item[0].data) rV = tuple(item[1].data) cV = tuple(item[2].data) rotMat = RotMat(rV) for i in range(numChildPerParent): nodeLabels.append(nodeLabel) nodeLabel += 1 arm = childParticleArm[i] armCurX = arm[0] * rotMat[0,0] + arm[1] * rotMat[0,1] + arm[2] * rotMat[0,2] armCurY = arm[0] * rotMat[1,0] + arm[1] * rotMat[1,1] + arm[2] * rotMat[1,2] armCurZ = arm[0] * rotMat[2,0] + arm[1] * rotMat[2,1] + arm[2] * rotMat[2,2] dispC0 = disp[0] + armCurX - arm[0] dispC1 = disp[1] + armCurY - arm[1] dispC2 = disp[2] + armCurZ - arm[2] cData = [] dispC = [] dispC.append(dispC0) dispC.append(dispC1) dispC.append(dispC2) cData0 = cV[0] + armCurX cData1 = cV[1] + armCurY cData2 = cV[2] + armCurZ cData.append(cData0) cData.append(cData1) cData.append(cData2) dispData.append(dispC) coordData.append(cData) else: for item in zip(uFieldSub.values, uRFieldSub.values): disp = tuple(item[0].data) rV = tuple(item[1].data) rotMat = RotMat(rV) for i in range(numChildPerParent): nodeLabels.append(nodeLabel) nodeLabel += 1 arm = childParticleArm[i] armCurX = arm[0] * rotMat[0,0] + arm[1] * rotMat[0,1] + arm[2] * rotMat[0,2] armCurY = arm[0] * rotMat[1,0] + arm[1] * rotMat[1,1] + arm[2] * rotMat[1,2] armCurZ = arm[0] * rotMat[2,0] + arm[1] * rotMat[2,1] + arm[2] * rotMat[2,2] dispC0 = disp[0] + armCurX - arm[0] dispC1 = disp[1] + armCurY - arm[1] dispC2 = disp[2] + armCurZ - arm[2] dispC = [] dispC.append(dispC0) dispC.append(dispC1) dispC.append(dispC2) dispData.append(dispC) # nodeLabel = nodeStartLabel for val in pradFieldSub.values: for i in range(numChildPerParent): Rad1 = (ChildParticleRad[i],) nodeLabel += 1 pradData.append(Rad1) # else: nodeLabel = nodeStartLabel if (iscField): for item in zip(uFieldSub.values, uRFieldSub.values, cFieldSub.values): disp = tuple(item[0].data) rV = tuple(item[1].data) cV = tuple(item[2].data) rotMat = RotMat(rV) for i in range(numChildPerParent): arm = childParticleArm[i] nodeLabels.append(nodeLabel) nodeLabel += 1 armCurX = arm[0] * rotMat[0,0] + arm[1] * rotMat[0,1] + arm[2] * rotMat[0,2] armCurY = arm[0] * rotMat[1,0] + arm[1] * rotMat[1,1] + arm[2] * rotMat[1,2] armCurZ = arm[0] * rotMat[2,0] + arm[1] * rotMat[2,1] + arm[2] * rotMat[2,2] cData = [] dispC = [] dispC0 = disp[0] + armCurX - arm[0] dispC1 = disp[1] + armCurY - arm[1] dispC2 = disp[2] + armCurZ - arm[2] dispC.append(dispC0) dispC.append(dispC1) dispC.append(dispC2) cData0 = cV[0] + armCurX cData1 = cV[1] + armCurY cData2 = cV[2] + armCurZ cData.append(cData0) cData.append(cData1) cData.append(cData2) dispData.append(dispC) coordData.append(cData) else: for item in zip(uFieldSub.values, uRFieldSub.values): disp = tuple(item[0].data) rV = tuple(item[1].data) rotMat = RotMat(rV) for i in range(numChildPerParent): arm = childParticleArm[i] nodeLabels.append(nodeLabel) nodeLabel += 1 armCurX = arm[0] * rotMat[0,0] + arm[1] * rotMat[0,1] + arm[2] * rotMat[0,2] armCurY = arm[0] * rotMat[1,0] + arm[1] * rotMat[1,1] + arm[2] * rotMat[1,2] armCurZ = arm[0] * rotMat[2,0] + arm[1] * rotMat[2,1] + arm[2] * rotMat[2,2] dispC = [] dispC0 = disp[0] + armCurX - arm[0] dispC1 = disp[1] + armCurY - arm[1] dispC2 = disp[2] + armCurZ - arm[2] dispC.append(dispC0) dispC.append(dispC1) dispC.append(dispC2) dispData.append(dispC) # nodeLabel = nodeStartLabel for val in pradFieldSub.values: for i in range(numChildPerParent): Rad1 = (ChildParticleRad[i],) nodeLabel += 1 pradData.append(Rad1) # uField.addData(position=NODAL, instance=inst, labels=nodeLabels, data=dispData) if (iscField): cField.addData(position=NODAL, instance=inst, labels=nodeLabels, data=coordData) pradField.addData(position=NODAL, instance=inst, labels=nodeLabels, data=pradData) #------------------------------------------------------------------------------- # Utility function to extract rotation matrix def RotMat(rV): rVMag2 = rV[0]*rV[0] + rV[1]*rV[1] + rV[2]*rV[2] rVMag = np.sqrt(rVMag2) SinTheta = math.sin(rVMag) CosTheta = math.cos(rVMag) if rVMag > 0.0: S1 = SinTheta/rVMag C1 = (1.0 - CosTheta)/(rVMag2) else: S1 = 0.0 C1 = 1.0 # SkwM = np.zeros((3,3)) #print(SkwM) SkwM[0][1] = -rV[2] SkwM[0][2] = rV[1] SkwM[1][0] = rV[2] SkwM[1][2] = -rV[0] SkwM[2][0] = -rV[1] SkwM[2][1] = rV[0] # rVrV = np.zeros((3,3)) rVrV[0][0] = rV[0]*rV[0] rVrV[0][1] = rV[0]*rV[1] rVrV[0][2] = rV[0]*rV[2] rVrV[1][0] = rV[1]*rV[0] rVrV[1][1] = rV[1]*rV[1] rVrV[1][2] = rV[1]*rV[2] rVrV[2][0] = rV[2]*rV[0] rVrV[2][1] = rV[2]*rV[1] rVrV[2][2] = rV[2]*rV[2] # RotM = np.zeros((3,3)) RotM[0][0] = CosTheta + S1*SkwM[0][0] + C1*rVrV[0][0] RotM[0][1] = S1*SkwM[0][1] + C1*rVrV[0][1] RotM[0][2] = S1*SkwM[0][2] + C1*rVrV[0][2] RotM[1][0] = S1*SkwM[1][0] + C1*rVrV[1][0] RotM[1][1] = CosTheta + S1*SkwM[1][1] + C1*rVrV[1][1] RotM[1][2] = S1*SkwM[1][2] + C1*rVrV[1][2] RotM[2][0] = S1*SkwM[2][0] + C1*rVrV[2][0] RotM[2][1] = S1*SkwM[2][1] + C1*rVrV[2][1] RotM[2][2] = CosTheta + S1*SkwM[2][2] + C1*rVrV[2][2] return RotM #------------------------------------------------------------------------------- # Utility routine to extract user specified orientation of the cluster def ClstOri(xAxis, yAxis, zAxis, vec): one = 1.0 vxx = xAxis[0] vyx = xAxis[1] vzx = xAxis[2] anorm = one / sqrt ( vxx * vxx + vyx * vyx + vzx * vzx ) vxx = anorm * vxx vyx = anorm * vyx vzx = anorm * vzx # vxz = vyx * yAxis[2] - vzx * yAxis[1] vyz = vzx * yAxis[0] - vxx * yAxis[2] vzz = vxx * yAxis[1] - vyx * yAxis[0] anorm = one / sqrt ( vxz * vxz + vyz * vyz + vzz * vzz ) vxz = anorm * vxz vyz = anorm * vyz vzz = anorm * vzz # vxy = vyz * vzx - vzz * vyx vyy = vzz * vxx - vxz * vzx vzy = vxz * vyx - vyz * vxx anorm = one / sqrt ( vxy * vxy + vyy * vyy + vzy * vzy ) vxy = anorm * vxy vyy = anorm * vyy vzy = anorm * vzy # vec[0][0] = vxx vec[1][0] = vyx vec[2][0] = vzx vec[0][1] = vxy vec[1][1] = vyy vec[2][1] = vzy vec[0][2] = vxz vec[1][2] = vyz vec[2][2] = vzz return vec print("Complete - copying the Odb...at", time.time()) odb = odbAccess.openOdb(dst, readOnly = False) if odb.isReadOnly == True: print('Error: odb is read only, check for .lck file and delete') exit() #------------------------------------------------------------------------------- # main driver section discreteSectionsInOdb = [] discreteSections = {} discreteSectionAssignments = {} numDiscreteSections = 0 sections = odb.sections def printSectionAssignment(sa): print("Section assignment: ") print(" section: ", sa.sectionName) print(" region name: ", sa.region.name) rp = 0 allSecChildParticleList = [] allSecRp = [] allSecElemSetName = [] for secName,section in sections.items(): if type(section) == DiscreteSectionType: discreteSections[secName] = section a = odb.rootAssembly for piname,pinstance in a.instances.items(): for sa in pinstance.sectionAssignments: if sa.sectionName in discreteSections.keys(): discreteSectionAssignments[sa.sectionName] = sa discreteSectionsInOdb.append(sa) for dsecName,dsection in discreteSections.items(): radiusType = dsection.radiusType if radiusType == SymbolicConstant("UNIFORM"): rp = dsection.radius allSecRp.append(rp) numDiscreteSections = numDiscreteSections + 1 childParticleList = [] childParticleData = dsection.childParticleData szChildParticleData = int(len(childParticleData) / 4) print(" number child particles per parent: %d" % szChildParticleData) for id in range(0,szChildParticleData): idx = id*4 lines = [] lines.append(childParticleData[idx]) lines.append(childParticleData[idx+1]) lines.append(childParticleData[idx+2]) lines.append(childParticleData[idx+3]) childParticleList.append(lines) allSecChildParticleList.append(childParticleList) if dsecName in discreteSectionAssignments: printSectionAssignment(discreteSectionAssignments[dsecName]) print("num Discrete sections:",numDiscreteSections) for id in range(numDiscreteSections): Rp = float(allSecRp[id]) section = discreteSectionsInOdb[id] createChildParticles(odb, str(id), section, Rp, allSecChildParticleList[id]) for id in range(numDiscreteSections): Rp = float(allSecRp[id]) section = discreteSectionsInOdb[id] parentSection = discreteSectionsInOdb[id] childInstanceName = parentSection.sectionName + '_child-1' applyDisplacement(odb, childInstanceName, parentSection, childParticleRadiusSection[id], childParticleArmSection[id] ) print('Done ...at', time.time()) print('Closing the odb...') odb.save() odb.close() print('Total time taken: %d sec'%(time.time()-startTime))