Location: Human Stomach Scaffold @ b24992e08081 / LoadxiInvertedmesh.py

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.cellml.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)