- Author:
- rjag008 <rjag008@auckland.ac.nz>
- Date:
- 2018-06-07 14:22:45+12:00
- Desc:
- Updated Provenance workflow uri
- Permanent Source URI:
- https://models.physiomeproject.org/workspace/516/rawfile/b24992e080818ebc5c9094afd2c38ec2d100cf0e/LoadxiInvertedmesh.py
'''
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)
