diff --git a/scaffoldmaker/utils/eftfactory_bicubichermitelinear.py b/scaffoldmaker/utils/eftfactory_bicubichermitelinear.py index 624c69e2..0d915bbb 100644 --- a/scaffoldmaker/utils/eftfactory_bicubichermitelinear.py +++ b/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 diff --git a/scaffoldmaker/utils/tubemesh.py b/scaffoldmaker/utils/tubemesh.py index 77b1a659..b7d8f404 100644 --- a/scaffoldmaker/utils/tubemesh.py +++ b/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 diff --git a/scaffoldmaker/utils/zinc_utils.py b/scaffoldmaker/utils/zinc_utils.py index 6c6fe294..002547ed 100644 --- a/scaffoldmaker/utils/zinc_utils.py +++ b/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