''' Created on 5/06/2018 @author: rjag008 ''' from __future__ import print_function from opencmiss.zinc.context import Context import numpy as np from opencmiss.zinc.node import Node from opencmiss.zinc.element import Elementbasis, Element class Stomach(object): ''' Loads stomach mesh generated by Kumar.. the element number is opposite to that of xi1 The generated mesh does not support cross derivatives ''' def __init__(self,ex2file,circumferentialElements=8,axialElements=11,wallElements=3): ''' Ex2 file containing node, element and field information ''' self.ex2 = ex2file self.circumferentialElements = circumferentialElements self.axialElements = axialElements self.wallElements = wallElements def normalizeToUnitBoundingVolume(self,filename): ''' Normalize the mesh to lie within a unit volume cube and centered around the origin ''' ctx = Context('Load') region = ctx.getDefaultRegion() region.readFile(self.ex2) fieldModule = region.getFieldmodule() coordinatesField = fieldModule.findFieldByName('coordinates') nodeset = fieldModule.findNodesetByFieldDomainType(coordinatesField.DOMAIN_TYPE_NODES) nodeIterator = nodeset.createNodeiterator () nodes = dict() node = nodeIterator.next() while node.isValid(): nodes[int(node.getIdentifier())-1] = node node = nodeIterator.next() fieldCache = fieldModule.createFieldcache() coordinates = np.zeros((len(nodes),3)) for nd,node in nodes.items(): fieldCache.setNode(node) _,coord = coordinatesField.evaluateReal(fieldCache,3) coordinates[nd-1] = coord minc = np.min(coordinates,axis=0) maxc = np.max(coordinates,axis=0) centroid = np.mean(coordinates,axis=0) normCoord = (coordinates-centroid)/(maxc-minc) fieldModule.beginChange() for nd,node in nodes.items(): fieldCache.setNode(node) coordinatesField.assignReal(fieldCache,list(normCoord[nd-1])) fieldModule.endChange() smoothing = fieldModule.createFieldsmoothing() coordinatesField.smooth(smoothing) region.writeFile(filename) def loadFromJson(self,filename): import json with open(filename,'r') as ser: result = json.load(ser) self.circumferentialElements = result['CircumferentialElements'] self.axialElements = result['AxialElements'] self.wallElements = result['WallElements'] if result['crossDerivatives']: raise ValueError('Cross derivatives are not suppported!') coordinates = np.array(result['coordinates']) fibres = np.array(result['fibres']) self.context = Context('HostMesh') region = self.context.getDefaultRegion() nodes,elems = self.generateTube(region, self.circumferentialElements, self.axialElements, self.wallElements) fieldModule = region.getFieldmodule() coordinatesField = fieldModule.findFieldByName('coordinates').castFiniteElement() fibresField = fieldModule.findFieldByName('fibres').castFiniteElement() fieldCache = fieldModule.createFieldcache() cubicHermiteNodeValuesPerDim = [Node.VALUE_LABEL_VALUE,Node.VALUE_LABEL_D_DS1,Node.VALUE_LABEL_D_DS2,\ Node.VALUE_LABEL_D_DS3] fieldModule.beginChange() for nd,node in nodes.items(): fieldCache.setNode(node) cv = coordinates[nd-1] fv = fibres[nd-1] for ct,nv in enumerate(cubicHermiteNodeValuesPerDim): coordinatesField.setNodeParameters(fieldCache,-1,nv,1,list(cv[:,ct])) fibresField.assignReal(fieldCache,list(fv)) fieldModule.endChange() region.writeFile('test.ex2') def generateJSON(self,filename): ctx = Context('Load') region = ctx.getDefaultRegion() region.readFile(self.ex2) fieldModule = region.getFieldmodule() coordinatesField = fieldModule.findFieldByName('coordinates').castFiniteElement() fibresField = fieldModule.findFieldByName('fibres').castFiniteElement() nodeset = fieldModule.findNodesetByFieldDomainType(coordinatesField.DOMAIN_TYPE_NODES) nodeIterator = nodeset.createNodeiterator () nodes = dict() node = nodeIterator.next() while node.isValid(): nodes[int(node.getIdentifier())-1] = node node = nodeIterator.next() fieldCache = fieldModule.createFieldcache() coordinates = np.zeros((len(nodes),3,4)) fibres = np.zeros((len(nodes),3)) cubicHermiteNodeValuesPerDim = [Node.VALUE_LABEL_VALUE,Node.VALUE_LABEL_D_DS1,Node.VALUE_LABEL_D_DS2,\ Node.VALUE_LABEL_D_DS3] for nd,node in nodes.items(): fieldCache.setNode(node) for cn in [1,2,3]: for ct,nv in enumerate(cubicHermiteNodeValuesPerDim): _,v = coordinatesField.getNodeParameters(fieldCache,cn,nv,1,1) coordinates[nd,cn-1,ct] = v _,fibre = fibresField.evaluateReal(fieldCache,3) fibres[nd-1] = fibre result = dict() result['CircumferentialElements'] = self.circumferentialElements result['AxialElements'] = self.axialElements result['WallElements'] = self.wallElements result['crossDerivatives'] = False result['coordinates'] = coordinates.tolist() result['fibres'] = fibres.tolist() import json with open(filename,'w') as ser: json.dump(result,ser) def generateTube(self,region,circumferentialElements,axialElements,wallElements,wallThickness=1): fieldModule = region.getFieldmodule() fieldModule.beginChange() coordinates = fieldModule.createFieldFiniteElement(3) coordinates.setName('coordinates') coordinates.setManaged(True) coordinates.setTypeCoordinate(True) coordinates.setCoordinateSystemType(coordinates.COORDINATE_SYSTEM_TYPE_RECTANGULAR_CARTESIAN) coordinates.setComponentName(1, 'x') coordinates.setComponentName(2, 'y') coordinates.setComponentName(3, 'z') fibres = fieldModule.createFieldFiniteElement(3) fibres.setName('fibres') fibres.setComponentName(1, 'fibre angle') fibres.setComponentName(2, 'imbrication angle') fibres.setComponentName(3, 'sheet angle') nodeset = fieldModule.findNodesetByFieldDomainType(coordinates.DOMAIN_TYPE_NODES) nodetemplate = nodeset.createNodetemplate() nodetemplate.defineField(coordinates) nodetemplate.defineField(fibres) for field in [coordinates]: nodetemplate.setValueNumberOfVersions(field, -1, Node.VALUE_LABEL_VALUE, 1) nodetemplate.setValueNumberOfVersions(field, -1, Node.VALUE_LABEL_D_DS1, 1) nodetemplate.setValueNumberOfVersions(field, -1, Node.VALUE_LABEL_D_DS2, 1) nodetemplate.setValueNumberOfVersions(field, -1, Node.VALUE_LABEL_D_DS3, 1) mesh = fieldModule.findMeshByDimension(3) tricubicHermiteBasis = fieldModule.createElementbasis(3, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE) eft = mesh.createElementfieldtemplate(tricubicHermiteBasis) #Setup the template such that the cross derivative terms are zero for n in range(8): eft.setFunctionNumberOfTerms(n*8 + 4, 0) eft.setFunctionNumberOfTerms(n*8 + 6, 0) eft.setFunctionNumberOfTerms(n*8 + 7, 0) eft.setFunctionNumberOfTerms(n*8 + 8, 0) elementtemplate = mesh.createElementtemplate() elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) result = elementtemplate.defineField(coordinates, -1, eft) result = elementtemplate.defineField(fibres, -1, eft) fieldCache = fieldModule.createFieldcache() nodes = dict() # create nodes radiansPerElementAround = 2.0*np.pi/circumferentialElements wallThicknessPerElement = wallThickness/wallElements x = [ 0.0, 0.0, 0.0 ] dx_ds1 = [ 0.0, 0.0, 0.0 ] dx_ds2 = [ 0.0, 0.0, 1.0 / axialElements ] dx_ds3 = [ 0.0, 0.0, 0.0 ] numberOfWallNodes = wallElements + 1 numberOfCircumfrentialNodes = circumferentialElements numberOfLengthNodes = axialElements + 1 for wallNodeIdx in range(1,numberOfWallNodes+1): radius = 0.5 + wallThickness*((wallNodeIdx-1)/(numberOfWallNodes - 1.0)) for lengthNodeIdx in range(1,numberOfLengthNodes+1): x[2] = float(lengthNodeIdx-1)/ axialElements for circumfrentialNodeIdx in range(1,numberOfCircumfrentialNodes+1): nodeNumber = circumfrentialNodeIdx + (lengthNodeIdx-1)*numberOfCircumfrentialNodes + (wallNodeIdx-1)*numberOfCircumfrentialNodes*numberOfLengthNodes radiansAround = circumfrentialNodeIdx*radiansPerElementAround cosRadiansAround = np.cos(radiansAround) sinRadiansAround = np.sin(radiansAround) x[0] = radius*cosRadiansAround x[1] = radius*sinRadiansAround dx_ds1[0] = radiansPerElementAround*radius*-sinRadiansAround dx_ds1[1] = radiansPerElementAround*radius*cosRadiansAround dx_ds3[0] = wallThicknessPerElement*cosRadiansAround dx_ds3[1] = wallThicknessPerElement*sinRadiansAround node = nodeset.createNode(nodeNumber, nodetemplate) nodes[nodeNumber] = node fieldCache.setNode(node) for coord in [coordinates]: coord.setNodeParameters(fieldCache, -1, Node.VALUE_LABEL_VALUE, 1, x) coord.setNodeParameters(fieldCache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1) coord.setNodeParameters(fieldCache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2) coord.setNodeParameters(fieldCache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) # create elements elementNumber = 0 elems = dict() for wallElementIdx in range(1,wallElements+1): for lengthElementIdx in range(1,axialElements+1): for circumfrentialElementIdx in range(1,circumferentialElements+1): elementNumber = elementNumber + 1 localNode1 = circumfrentialElementIdx + (lengthElementIdx - 1)*circumferentialElements + \ (wallElementIdx-1)*numberOfCircumfrentialNodes*numberOfLengthNodes if circumfrentialElementIdx == circumferentialElements: localNode2 = 1 + (lengthElementIdx-1)*numberOfCircumfrentialNodes + \ (wallElementIdx-1)*numberOfCircumfrentialNodes*numberOfLengthNodes else: localNode2 = localNode1 + 1 localNode3 = localNode1 + numberOfCircumfrentialNodes localNode4 = localNode2 + numberOfCircumfrentialNodes localNode5 = localNode1 + numberOfCircumfrentialNodes*numberOfLengthNodes localNode6 = localNode2 + numberOfCircumfrentialNodes*numberOfLengthNodes localNode7 = localNode3 + numberOfCircumfrentialNodes*numberOfLengthNodes localNode8 = localNode4 + numberOfCircumfrentialNodes*numberOfLengthNodes localNodes = [localNode1,localNode2,localNode3,localNode4,localNode5,localNode6,localNode7,localNode8] element = mesh.createElement(elementNumber, elementtemplate) result = element.setNodesByIdentifier(eft, localNodes) elems[elementNumber] = element fieldModule.defineAllFaces() fieldModule.endChange() return nodes,elems def getInitialValues(self,fieldNames,circumferentialElements,axialElements,wallElements,refineAtLength,refineAtTheta): ''' Determine initial values of fields in fieldNames from source mesh based on material coordinates of target mesh nodes fieldNames is a dictionary with fieldName,numberOfComponents {'coordinates':3,'label',:1} ''' numberOfCircumfrentialNodes = circumferentialElements numberOfLengthNodes = axialElements + 1 # create node element map ev = 0 nodes = dict() linterval = axialElements-len(refineAtLength) if linterval > 8: lengthElementLocations = np.linspace(0.0, 0.99999, linterval+1) elif linterval < 4: raise ValueError('At least four elements are required along the axis') else: eloc= np.linspace(0.0, 1.0, 11) le = np.linspace(eloc[1], eloc[-2], linterval-1) le1 = [0.0] le1.extend(le.tolist()) le1.append(0.9999) lengthElementLocations = np.array(le1) cinterval = circumferentialElements-len(refineAtTheta) if cinterval < 4: raise ValueError('At least four elements are required along the circumference') circumferentialElementLocations = np.linspace(0.0, 0.99999, cinterval+1) ctr = 0 for elem,xi in refineAtTheta.items(): nxi = xi/float(cinterval)+circumferentialElementLocations[elem-1] circumferentialElementLocations = np.insert(circumferentialElementLocations,elem+ctr,nxi) ctr +=1 ctr = 0 for elem,xi in refineAtLength.items(): nxi = xi/float(linterval)+lengthElementLocations[elem-1] lengthElementLocations = np.insert(lengthElementLocations,elem+ctr,nxi) ctr +=1 wallElementLocations = np.linspace(0.0, 0.99999, wallElements+1) #Compute the material coordinates of nodes for wallElementIdx in range(1,wallElements+1): for lengthElementIdx in range(1,axialElements+1): for circumfrentialElementIdx in range(1,circumferentialElements+1): #Element coord with respect to elements along each dim localNode1 = circumfrentialElementIdx + (lengthElementIdx - 1)*circumferentialElements + \ (wallElementIdx-1)*numberOfCircumfrentialNodes*numberOfLengthNodes if circumfrentialElementIdx == circumferentialElements: localNode2 = 1 + (lengthElementIdx-1)*numberOfCircumfrentialNodes + \ (wallElementIdx-1)*numberOfCircumfrentialNodes*numberOfLengthNodes else: localNode2 = localNode1 + 1 localNode3 = localNode1 + numberOfCircumfrentialNodes localNode4 = localNode2 + numberOfCircumfrentialNodes localNode5 = localNode1 + numberOfCircumfrentialNodes*numberOfLengthNodes localNode6 = localNode2 + numberOfCircumfrentialNodes*numberOfLengthNodes localNode7 = localNode3 + numberOfCircumfrentialNodes*numberOfLengthNodes localNode8 = localNode4 + numberOfCircumfrentialNodes*numberOfLengthNodes l1xi = [circumferentialElementLocations[circumfrentialElementIdx-1],lengthElementLocations[lengthElementIdx-1],wallElementLocations[wallElementIdx-1]] l2xi = [circumferentialElementLocations[circumfrentialElementIdx],lengthElementLocations[lengthElementIdx-1],wallElementLocations[wallElementIdx-1]] l3xi = [circumferentialElementLocations[circumfrentialElementIdx-1],lengthElementLocations[lengthElementIdx],wallElementLocations[wallElementIdx-1]] l4xi = [circumferentialElementLocations[circumfrentialElementIdx],lengthElementLocations[lengthElementIdx],wallElementLocations[wallElementIdx-1]] l5xi = [circumferentialElementLocations[circumfrentialElementIdx-1],lengthElementLocations[lengthElementIdx-1],wallElementLocations[wallElementIdx]] l6xi = [circumferentialElementLocations[circumfrentialElementIdx],lengthElementLocations[lengthElementIdx-1],wallElementLocations[wallElementIdx]] l7xi = [circumferentialElementLocations[circumfrentialElementIdx-1],lengthElementLocations[lengthElementIdx],wallElementLocations[wallElementIdx]] l8xi = [circumferentialElementLocations[circumfrentialElementIdx],lengthElementLocations[lengthElementIdx],wallElementLocations[wallElementIdx]] ev +=1 nds = [localNode1,localNode2,localNode3,localNode4,localNode5,localNode6,localNode7,localNode8] xis = [l1xi,l2xi,l3xi,l4xi,l5xi,l6xi,l7xi,l8xi] for i,nd in enumerate(nds): if not nd in nodes: nodes[nd] = xis[i] ctx = Context('Load') region = ctx.getDefaultRegion() region.readFile(self.ex2) fieldModule = region.getFieldmodule() #Coordinates is a cubic hermite Field fieldHandles = dict() for field in fieldNames.keys(): hField = fieldModule.findFieldByName(field) hFielddx1 = fieldModule.createFieldDerivative(hField,1) hFielddx2 = fieldModule.createFieldDerivative(hField,2) hFielddx3 = fieldModule.createFieldDerivative(hField,3) fieldHandles[field]=[hField,hFielddx1,hFielddx2,hFielddx3] #Get the elements mesh = fieldModule.findMeshByDimension(3) ei = mesh.createElementiterator() elemDict = dict() elem = ei.next() while elem.isValid(): elemDict[int(elem.getIdentifier())] = elem elem = ei.next() def getLocationValuesAndDerivatives(fields): _,v = hfl[0].evaluateReal(fieldCache,fieldNames[fn]) _,dx1 = fields[1].evaluateReal(fieldCache,3)#d/ds1 _,dx2 = fields[2].evaluateReal(fieldCache,3)#d/ds2 _,dx3 = fields[3].evaluateReal(fieldCache,3)#d/ds3 val = [v,dx1,dx2,dx3] return val fieldCache = fieldModule.createFieldcache() for nd,val in nodes.items(): #cval = np.clip(val,0,0.9999) cval = val xpos = int(cval[0]*self.circumferentialElements) xi1 = val[0]*self.circumferentialElements - xpos ypos = int(cval[1]*self.axialElements) xi2 = val[1]*self.axialElements - ypos zpos = int(cval[2]*self.wallElements) xi3 = val[2]*self.wallElements - zpos eid = xpos + ypos*self.circumferentialElements + zpos*(self.circumferentialElements*self.axialElements) + 1 fv = dict() #xi1 and xi2 are inverted in kumar's mesh for fn,hfl in fieldHandles.items(): fieldCache.setMeshLocation(elemDict[eid],[xi2,xi1,xi3]) #_,v = hfl[0].evaluateReal(fieldCache,fieldNames[fn]) #fv[fn] = v fv[fn] = getLocationValuesAndDerivatives(hfl) nodes[nd] = fv return nodes def createMesh(self,filename="test.ex2",circumferentialElements=8,axialElements=11,wallElements=3,refineAtLength={},refineAtTheta={}): ''' filename - name of output file circumferentialElements - number of elements along the circumference axialElements - number of elements along the axis wallElements - number of elements along the wall refineAtLength - dict of axial length that need to be refined at an xi {1:0.25}, will use coordinates from 0-0.25 for 1 and 0.25-0.75 for 2 Additional elements along the axis will be added for each entry refineAtTheta - list of circumferential elements that need to be refined ''' ctx = Context('Load') region = ctx.getDefaultRegion() nodes,_ = self.generateTube(region, circumferentialElements+len(refineAtTheta), axialElements+len(refineAtLength), wallElements) nvals = self.getInitialValues({'coordinates':3}, circumferentialElements+len(refineAtTheta), axialElements+len(refineAtLength), \ wallElements,refineAtLength,refineAtTheta) fieldModule = region.getFieldmodule() fieldCache = fieldModule.createFieldcache() coordinatesField = fieldModule.findFieldByName('coordinates').castFiniteElement() fieldModule.beginChange() for nd,node in nodes.items(): cv = nvals[nd]['coordinates'] fieldCache.setNode(node) #coordinatesField.assignReal(fieldCache,cv) coordinatesField.setNodeParameters(fieldCache, -1, Node.VALUE_LABEL_VALUE, 1, cv[0]) #coordinatesField.setNodeParameters(fieldCache, -1, Node.VALUE_LABEL_D_DS1, 1, cv[1]) #coordinatesField.setNodeParameters(fieldCache, -1, Node.VALUE_LABEL_D_DS2, 1, cv[2]) #coordinatesField.setNodeParameters(fieldCache, -1, Node.VALUE_LABEL_D_DS3, 1, cv[3]) print("print Statement at LoadxiInverteredmesh 254: Using smooth until derivatives can be computed explicitly") smoothing = fieldModule.createFieldsmoothing() coordinatesField.smooth(smoothing) fieldModule.endChange() region.writeFile(filename) if __name__ == "__main__": obj = Stomach('finalstomach.ex2') obj.normalizeToUnitBoundingVolume('normalizedStomach.ex2') refineLatitute={6:0.5,1:0.25} #Refine elements at the current 6th and 1st axial number at xi 0.5 and 0.25 refineLongitude={1:0.5} #Refine elements at the current 1st circumferential number at xi 0.5 obj.createMesh(refineAtLength=refineLatitute,refineAtTheta=refineLongitude)