Merge pull request #46 from mlin865/textureCoord

Texture coordinates

scaffoldmaker/utils/eftfactory_bicubichermitelinear.py

@@ -206,3 +206,23 @@ def createEftSplitXi1RightOut(self):
         remapEftNodeValueLabel(eft, [ 5, 7 ], self._d_ds2, [ (self._d_ds1, [1]), (self._d_ds2, []) ])
         assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftSplitXi1RightOut:  Failed to validate eft'
         return eft
+
+    def createEftOpenTube(self):
+        '''
+        Create a basic bicubic hermite linear element template for elements
+        along boundary where a tube is opened on xi1 = 1 for a flat preparation.
+        Could eventually have 6 variants. Retain node numbering with two versions
+        for boundary nodes.
+        :return: Element field template
+        '''
+        eft = self.createEftBasic()
+        for n in [ 1, 3, 5, 7 ]:
+            ln = n + 1
+            eft.setTermNodeParameter(n*4 + 1, 1, ln, Node.VALUE_LABEL_VALUE, 2)
+            eft.setTermNodeParameter(n*4 + 2, 1, ln, Node.VALUE_LABEL_D_DS1, 2)
+            eft.setTermNodeParameter(n*4 + 3, 1, ln, Node.VALUE_LABEL_D_DS2, 2)
+            if self._useCrossDerivatives:
+                eft.setTermNodeParameter(n*4 + 4, 1, ln, Node.VALUE_LABEL_D2_DS1DS2, 2)
+
+        assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftOpenTube:  Failed to validate eft'
+        return eft

scaffoldmaker/utils/tubemesh.py

@@ -25,7 +25,7 @@ def generatetubemesh(region,
     segmentAxis, segmentLength,
     useCrossDerivatives,
     useCubicHermiteThroughWall, # or Zinc Elementbasis.FUNCTION_TYPE_LINEAR_LAGRANGE etc.
-    nextNodeIdentifier = 1, nextElementIdentifier = 1
+    firstNodeIdentifier = 1, firstElementIdentifier = 1
     ):
     '''
     Generates a 3-D tubular mesh with variable numbers of elements
@@ -45,7 +45,7 @@ def generatetubemesh(region,
     :param segmentAxis: axis of segment profile
     :param segmentLength: length of segment profile
     :param useCubicHermiteThroughWall: use linear when false
-    :return: annotationGroups, nextNodeIdentifier, nextElementIdentifier
+    :return: annotationGroups, nodeIdentifier, elementIdentifier
     '''
     zero  = [0.0, 0.0, 0.0]
     annotationGroups = []
@@ -71,7 +71,7 @@ def generatetubemesh(region,
 
     # Step through central line and rotate central line axes to align tangent 
     # to tangent from previous frame
-    for n in range(1, elementsCountAlong+1):
+    for n in range(1, elementsCountAlong + 1):
         unitTangent = normalise(sd1[n])
         cp = crossproduct3(prevUnitTangent, unitTangent)
         if magnitude(cp)> 0.0:
@@ -122,7 +122,7 @@ def generatetubemesh(region,
     result = elementtemplate.defineField(coordinates, -1, eft)
 
     # create nodes
-    nodeIdentifier = nextNodeIdentifier
+    nodeIdentifier = firstNodeIdentifier
     x = [ 0.0, 0.0, 0.0 ]
     dx_ds1 = [ 0.0, 0.0, 0.0 ]
     dx_ds2 = [ 0.0, 0.0, 0.0 ]
@@ -138,7 +138,7 @@ def generatetubemesh(region,
 
     # Map each face along segment profile to central line
     for nSegment in range(segmentCountAlong):
-        for nAlongSegment in range(elementsCountAlongSegment+1):
+        for nAlongSegment in range(elementsCountAlongSegment + 1):
             n2 = nSegment*elementsCountAlongSegment + nAlongSegment
             if nSegment == 0 or (nSegment > 0 and nAlongSegment > 0):
                 # Rotate to align segment axis with tangent of central line
@@ -229,20 +229,96 @@ def generatetubemesh(region,
         # nodeIdentifier = nodeIdentifier + 1
 
     # create elements
-    elementIdentifier = nextElementIdentifier
+    elementIdentifier = firstElementIdentifier
     now = (elementsCountAlong + 1)*elementsCountAround
     for e3 in range(elementsCountThroughWall):
         for e2 in range(elementsCountAlong):
             for e1 in range(elementsCountAround):
                 element = mesh.createElement(elementIdentifier, elementtemplate)
                 bni11 = e3*now + e2*elementsCountAround + e1 + 1
-                bni12 = e3*now + e2*elementsCountAround + (e1 + 1)%elementsCountAround + 1
+                bni12 = e3*now + e2*elementsCountAround + (e1 + 1) % elementsCountAround + 1
                 bni21 = e3*now + (e2 + 1)*elementsCountAround + e1 + 1
-                bni22 = e3*now + (e2 + 1)*elementsCountAround + (e1 + 1)%elementsCountAround + 1
+                bni22 = e3*now + (e2 + 1)*elementsCountAround + (e1 + 1) % elementsCountAround + 1
                 nodeIdentifiers = [ bni11, bni12, bni21, bni22, bni11 + now, bni12 + now, bni21 + now, bni22 + now ]
                 result = element.setNodesByIdentifier(eft, nodeIdentifiers)
                 elementIdentifier = elementIdentifier + 1
 
+    # Define texture coordinates field
+    textureCoordinates = getOrCreateTextureCoordinateField(fm)
+    textureNodetemplate1 = nodes.createNodetemplate()
+    textureNodetemplate1.defineField(textureCoordinates)
+    textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 1)
+    textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 1)
+    textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 1)
+    if useCrossDerivatives:
+        textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1)
+
+    textureNodetemplate2 = nodes.createNodetemplate()
+    textureNodetemplate2.defineField(textureCoordinates)
+    textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 2)
+    textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 2)
+    textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 2)
+    if useCrossDerivatives:
+        textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 2)
+
+    bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives)
+    eftTexture1 = bicubichermitelinear.createEftBasic()
+
+    elementtemplate1 = mesh.createElementtemplate()
+    elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE)
+    elementtemplate1.defineField(textureCoordinates, -1, eftTexture1)
+
+    eftTexture2 = bicubichermitelinear.createEftOpenTube()
+    elementtemplate2 = mesh.createElementtemplate()
+    elementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE)
+    elementtemplate2.defineField(textureCoordinates, -1, eftTexture2)
+
+    # Calculate texture coordinates and derivatives
+    d1 = [1.0 / elementsCountAround, 0.0, 0.0]
+    d2 = [0.0, 1.0 / elementsCountAlong, 0.0]
+
+    nodeIdentifier = firstNodeIdentifier
+    for n3 in range(elementsCountThroughWall + 1):
+        for n2 in range(elementsCountAlong + 1):
+            for n1 in range(elementsCountAround):
+                u = [ 1.0 / elementsCountAround * n1,
+                    1.0 / elementsCountAlong * n2,
+                    1.0 / elementsCountThroughWall * n3]
+                node = nodes.findNodeByIdentifier(nodeIdentifier)
+                node.merge(textureNodetemplate2 if n1 == 0 else textureNodetemplate1)
+                cache.setNode(node)
+                textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, u)
+                textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1)
+                textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2)
+                if useCrossDerivatives:
+                    textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero)
+                if n1 == 0:
+                    u = [ 1.0, 1.0 / elementsCountAlong * n2, 1.0 / elementsCountThroughWall * n3]
+                    textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, u)
+                    textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1)
+                    textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2)
+                    if useCrossDerivatives:
+                        textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero)
+                nodeIdentifier = nodeIdentifier + 1
+
+    # define texture coordinates field over elements
+    elementIdentifier = firstElementIdentifier
+    now = (elementsCountAlong + 1)*elementsCountAround
+
+    for e3 in range(elementsCountThroughWall):
+        for e2 in range(elementsCountAlong):
+            for e1 in range(elementsCountAround):
+                onOpening = e1 > elementsCountAround - 2
+                element = mesh.findElementByIdentifier(elementIdentifier)
+                element.merge(elementtemplate2 if onOpening else elementtemplate1)
+                bni11 = e3*now + e2*elementsCountAround + e1 + 1
+                bni12 = e3*now + e2*elementsCountAround + (e1 + 1) % elementsCountAround + 1
+                bni21 = e3*now + (e2 + 1)*elementsCountAround + e1 + 1
+                bni22 = e3*now + (e2 + 1)*elementsCountAround + (e1 + 1) % elementsCountAround + 1
+                nodeIdentifiers = [ bni11, bni12, bni21, bni22, bni11 + now, bni12 + now, bni21 + now, bni22 + now ]
+                element.setNodesByIdentifier(eftTexture2 if onOpening else eftTexture1, nodeIdentifiers)
+                elementIdentifier = elementIdentifier + 1
+
     fm.endChange()
 
     return annotationGroups, nodeIdentifier, elementIdentifier
@@ -267,8 +343,8 @@ def getOuterCoordinatesAndCurvatureFromInner(xInner, d1Inner, d3Inner, wallThick
         for n1 in range(elementsCountAround):
             n = n2*elementsCountAround + n1
             x = [xInner[n][i] + d3Inner[n][i]*wallThickness for i in range(3)]
-            prevIdx = n-1 if (n1 != 0) else (n2+1)*elementsCountAround - 1
-            nextIdx = n+1 if (n1 < elementsCountAround-1) else n2*elementsCountAround
+            prevIdx = n - 1 if (n1 != 0) else (n2 + 1)*elementsCountAround - 1
+            nextIdx = n + 1 if (n1 < elementsCountAround - 1) else n2*elementsCountAround
             norm = d3Inner[n]
             curvatureAround = 0.5*(
                 getCubicHermiteCurvature(xInner[prevIdx], d1Inner[prevIdx], xInner[n], d1Inner[n], norm, 1.0) +
@@ -300,7 +376,7 @@ def interpolatefromInnerAndOuter( xInner, xOuter, thickness, xi3, curvatureInner
     dx_ds2List = []
     dx_ds3List =[]
 
-    for n2 in range(elementsCountAlong+1):
+    for n2 in range(elementsCountAlong + 1):
         for n1 in range(elementsCountAround):
             n = n2*elementsCountAround + n1
             norm = d3InnerUnit[n]
@@ -316,16 +392,16 @@ def interpolatefromInnerAndOuter( xInner, xOuter, thickness, xi3, curvatureInner
             dx_ds1List.append(dx_ds1)
             # dx_ds2
             if n2 > 0 and n2 < elementsCountAlong:
-                prevIdx = (n2-1)*elementsCountAround + n1
-                nextIdx = (n2+1)*elementsCountAround + n1
+                prevIdx = (n2 - 1)*elementsCountAround + n1
+                nextIdx = (n2 + 1)*elementsCountAround + n1
                 curvatureAround = 0.5*(
-                    getCubicHermiteCurvature(xInner[prevIdx], d2Inner[prevIdx], xInner[n], d2Inner[n], norm, 1.0)+
+                    getCubicHermiteCurvature(xInner[prevIdx], d2Inner[prevIdx], xInner[n], d2Inner[n], norm, 1.0) +
                     getCubicHermiteCurvature(xInner[n], d2Inner[n], xInner[nextIdx], d2Inner[nextIdx], norm, 0.0))
             elif n2 == 0:
-                nextIdx = (n2+1)*elementsCountAround + n1
+                nextIdx = (n2 + 1)*elementsCountAround + n1
                 curvatureAround = getCubicHermiteCurvature(xInner[n], d2Inner[n], xInner[nextIdx], d2Inner[nextIdx], norm, 0.0)
             else:
-                prevIdx = (n2-1)*elementsCountAround + n1
+                prevIdx = (n2 - 1)*elementsCountAround + n1
                 curvatureAround = getCubicHermiteCurvature(xInner[prevIdx], d2Inner[prevIdx], xInner[n], d2Inner[n], norm, 1.0)
 
             factor = 1.0 - curvatureAround*thickness*xi3

scaffoldmaker/utils/zinc_utils.py

@@ -38,6 +38,34 @@ def getOrCreateCoordinateField(fieldmodule, name='coordinates', componentsCount=
     fieldmodule.endChange()
     return coordinates
 
+def getOrCreateTextureCoordinateField(fieldmodule, name='textureCoordinates', componentsCount=3):
+    '''
+    Finds or creates a rectangular cartesian texture coordinate field.
+    New field has component names, 'u', 'v', 'w'.
+    Raises exception if existing field of name is not finite element type or has incorrect attributes.
+    :param fieldmodule:  Zinc fieldmodule to find or create field in.
+    :param name:  Name of field to find or create.
+    :param componentsCount: Number of components / dimension of field.
+    '''
+    assert (componentsCount > 0) and (componentsCount <= 3), 'getOrCreateCoordinateField.  Dimensions must be from 1 to 3'
+    coordinates = fieldmodule.findFieldByName(name)
+    if coordinates.isValid():
+        coordinates = coordinates.castFiniteElement()
+        assert coordinates.isValid(), 'getOrCreateCoordinateField.  Existing field \'' + name + '\' is not finite element type'
+        assert coordinates.getNumberOfComponents() == componentsCount, 'getOrCreateCoordinateField.  Existing field \'' + name + '\' does not have ' + str(componentsCount) + ' components'
+        assert coordinates.getCoordinateSystemType() == Field.COORDINATE_SYSTEM_TYPE_RECTANGULAR_CARTESIAN, 'getOrCreateCoordinateField.  Existing field \'' + name + '\' is not rectangular Cartesian'
+        return coordinates
+    fieldmodule.beginChange()
+    coordinates = fieldmodule.createFieldFiniteElement(componentsCount)
+    coordinates.setName(name)
+    coordinates.setManaged(True)
+    coordinates.setTypeCoordinate(True)
+    coordinates.setCoordinateSystemType(Field.COORDINATE_SYSTEM_TYPE_RECTANGULAR_CARTESIAN)
+    for c in range(componentsCount):
+        coordinates.setComponentName(c + 1, ['u', 'v', 'w'][c])
+    fieldmodule.endChange()
+    return coordinates
+
 def getOrCreateElementXiField(fieldmodule, name='element_xi', mesh=None):
     '''
     Finds or creates a stored mesh location field for storing locations in the