diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py index 74ba22c5..3a5c9886 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py @@ -50,15 +50,16 @@ def getDefaultOptions(cls, parameterSetName="Default"): options["Head length"] = 2.2 options["Head width"] = 2.0 options["Neck length"] = 1.3 - options["Shoulder drop"] = 0.7 + options["Shoulder drop"] = 1.0 options["Shoulder width"] = 4.5 options["Arm lateral angle degrees"] = 10.0 options["Arm length"] = 7.5 options["Arm top diameter"] = 1.0 + options["Arm twist angle degrees"] = 0.0 options["Wrist thickness"] = 0.5 options["Wrist width"] = 0.7 - options["Hand length"] = 1.5 - options["Hand thickness"] = 0.3 + options["Hand length"] = 2.0 + options["Hand thickness"] = 0.2 options["Hand width"] = 1.0 options["Thorax length"] = 2.5 options["Abdomen length"] = 3.0 @@ -90,6 +91,7 @@ def getOrderedOptionNames(cls): "Arm lateral angle degrees", "Arm length", "Arm top diameter", + "Arm twist angle degrees", "Wrist thickness", "Wrist width", "Hand length", @@ -154,20 +156,15 @@ def checkOptions(cls, options): options[key] = 0.1 elif options[key] > 0.9: options[key] = 0.9 - for key in [ - "Arm lateral angle degrees" - ]: - if options[key] < -60.0: - options[key] = -60.0 - elif options[key] > 200.0: - options[key] = 200.0 - for key in [ - "Leg lateral angle degrees" - ]: - if options[key] < -20.0: - options[key] = -20.0 - elif options[key] > 60.0: - options[key] = 60.0 + for key, angleRange in { + "Arm lateral angle degrees": (-60.0, 200.0), + "Arm twist angle degrees": (-90.0, 90.0), + "Leg lateral angle degrees": (-20.0, 60.0) + }.items(): + if options[key] < angleRange[0]: + options[key] = angleRange[0] + elif options[key] > angleRange[1]: + options[key] = angleRange[1] return dependentChanges @classmethod @@ -189,6 +186,7 @@ def generateBaseMesh(cls, region, options): armAngleRadians = math.radians(options["Arm lateral angle degrees"]) armLength = options["Arm length"] armTopRadius = 0.5 * options["Arm top diameter"] + armTwistAngleRadians = math.radians(options["Arm twist angle degrees"]) halfWristThickness = 0.5 * options["Wrist thickness"] halfWristWidth = 0.5 * options["Wrist width"] handLength = options["Hand length"] @@ -396,22 +394,23 @@ def generateBaseMesh(cls, region, options): # initial shoulder rotation with arm is negligible, hence: shoulderRotationFactor = 1.0 - math.cos(0.5 * armAngleRadians) # assume shoulder drop is half shrug distance to get limiting shoulder angle for 180 degree arm rotation - shoulderLimitAngleRadians = math.asin(2.0 * shoulderDrop / halfShoulderWidth) + shoulderLimitAngleRadians = math.asin(1.5 * shoulderDrop / halfShoulderWidth) shoulderAngleRadians = shoulderRotationFactor * shoulderLimitAngleRadians armStartX = thoraxStartX + shoulderDrop - halfShoulderWidth * math.sin(shoulderAngleRadians) nonHandArmLength = armLength - handLength armScale = nonHandArmLength / (armToHandElementsCount - 2) # 2 == shoulder elements count d12_mag = (halfWristThickness - armTopRadius) / (armToHandElementsCount - 2) d13_mag = (halfWristWidth - armTopRadius) / (armToHandElementsCount - 2) - hd3 = [0.0, 0.0, halfHandWidth] - hid3 = mult(hd3, innerProportionDefault) for side in (left, right): armAngle = armAngleRadians if (side == left) else -armAngleRadians cosArmAngle = math.cos(armAngle) sinArmAngle = math.sin(armAngle) armStartY = (halfShoulderWidth if (side == left) else -halfShoulderWidth) * math.cos(shoulderAngleRadians) x = [armStartX, armStartY, 0.0] - d1 = [armScale * cosArmAngle, armScale * sinArmAngle, 0.0] + armDirn = [cosArmAngle, sinArmAngle, 0.0] + armSide = [-sinArmAngle, cosArmAngle, 0.0] + armFront = cross(armDirn, armSide) + d1 = mult(armDirn, armScale) # set leg versions 2 (left) and 3 (right) on leg junction node, and intermediate shoulder node sd1 = interpolateLagrangeHermiteDerivative(sx, x, d1, 0.0) nx, nd1 = sampleCubicHermiteCurvesSmooth([sx, x], [sd1, d1], 2, derivativeMagnitudeEnd=armScale)[0:2] @@ -454,11 +453,9 @@ def generateBaseMesh(cls, region, options): sid13 = mult(sd13, innerProportionDefault) innerCoordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, version, sid12) innerCoordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, version, sid13) - d12 = [-d12_mag * sinArmAngle, d12_mag * cosArmAngle, 0.0] - id12 = mult(d12, innerProportionDefault) - d13 = [0.0, 0.0, d13_mag] - id13 = mult(d13, innerProportionDefault) # main part of arm to wrist + elementTwistAngle = ((armTwistAngleRadians if (side == left) else -armTwistAngleRadians) / + (armToHandElementsCount - 3)) for i in range(armToHandElementsCount - 1): xi = i / (armToHandElementsCount - 2) node = nodes.findNodeByIdentifier(nodeIdentifier) @@ -466,10 +463,33 @@ def generateBaseMesh(cls, region, options): x = [armStartX + d1[0] * i, armStartY + d1[1] * i, d1[2] * i] halfThickness = xi * halfWristThickness + (1.0 - xi) * armTopRadius halfWidth = xi * halfWristWidth + (1.0 - xi) * armTopRadius - d2 = [-halfThickness * sinArmAngle, halfThickness * cosArmAngle, 0.0] - d3 = [0.0, 0.0, halfWidth] + if i == 0: + twistAngle = 0.0 + elif i == (armToHandElementsCount - 2): + twistAngle = armTwistAngleRadians if (side == left) else -armTwistAngleRadians + else: + twistAngle = -0.5 * elementTwistAngle + elementTwistAngle * i + if twistAngle == 0.0: + d2 = mult(armSide, halfThickness) + d3 = mult(armFront, halfWidth) + d12 = mult(armSide, d12_mag) + d13 = mult(armFront, d13_mag) + else: + cosTwistAngle = math.cos(twistAngle) + sinTwistAngle = math.sin(twistAngle) + d2 = sub(mult(armSide, halfThickness * cosTwistAngle), + mult(armFront, halfThickness * sinTwistAngle)) + d3 = add(mult(armFront, halfWidth * cosTwistAngle), + mult(armSide, halfWidth * sinTwistAngle)) + d12 = set_magnitude(d2, d12_mag) + d13 = set_magnitude(d3, d13_mag) + if i < (armToHandElementsCount - 2): + d12 = add(d12, set_magnitude(d3, -halfThickness * elementTwistAngle)) + d13 = add(d13, set_magnitude(d2, halfWidth * elementTwistAngle)) id2 = mult(d2, innerProportionDefault) id3 = mult(d3, innerProportionDefault) + id12 = mult(d12, innerProportionDefault) + id13 = mult(d13, innerProportionDefault) setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) nodeIdentifier += 1 @@ -479,8 +499,19 @@ def generateBaseMesh(cls, region, options): fieldcache.setNode(node) hx = [armStartX + armLength * cosArmAngle, armStartY + armLength * sinArmAngle, 0.0] hd1 = computeCubicHermiteEndDerivative(x, d1, hx, d1) - hd2 = set_magnitude(d2, halfHandThickness) + twistAngle = armTwistAngleRadians if (side == left) else -armTwistAngleRadians + if twistAngle == 0.0: + hd2 = set_magnitude(d2, halfHandThickness) + hd3 = [0.0, 0.0, halfHandWidth] + else: + cosTwistAngle = math.cos(twistAngle) + sinTwistAngle = math.sin(twistAngle) + hd2 = sub(mult(armSide, halfHandThickness * cosTwistAngle), + mult(armFront, halfHandThickness * sinTwistAngle)) + hd3 = add(mult(armFront, halfHandWidth * cosTwistAngle), + mult(armSide, halfHandWidth * sinTwistAngle)) hid2 = mult(hd2, innerProportionDefault) + hid3 = mult(hd3, innerProportionDefault) setNodeFieldParameters(coordinates, fieldcache, hx, hd1, hd2, hd3) setNodeFieldParameters(innerCoordinates, fieldcache, hx, hd1, hid2, hid3) nodeIdentifier += 1 diff --git a/tests/test_wholebody2.py b/tests/test_wholebody2.py index 69294fed..f59f448d 100644 --- a/tests/test_wholebody2.py +++ b/tests/test_wholebody2.py @@ -69,8 +69,8 @@ def test_wholebody2_core(self): self.assertTrue(coordinates.isValid()) minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) tol = 1.0E-4 - assertAlmostEqualList(self, minimums, [0.0, -3.7000751482231564, -1.25], tol) - assertAlmostEqualList(self, maximums, [20.48332368641937, 3.7000751482231564, 2.15], tol) + assertAlmostEqualList(self, minimums, [0.0, -3.650833433150559, -1.25], tol) + assertAlmostEqualList(self, maximums, [20.48332368641937, 3.650833433150559, 2.15], tol) with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) @@ -88,17 +88,17 @@ def test_wholebody2_core(self): result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 98.86108394657288, delta=tol) - self.assertAlmostEqual(surfaceArea, 229.91482811362502, delta=tol) + self.assertAlmostEqual(volume, 98.41184838749453, delta=tol) + self.assertAlmostEqual(surfaceArea, 228.9017722286146, delta=tol) # check some annotation groups: expectedSizes3d = { "abdominal cavity": (40, 10.074522362520469), - "core": (428, 45.718809433770865), + "core": (428, 45.535080392793354), "head": (64, 6.909618374858558), - "thoracic cavity": (40, 7.135491643165788), - "shell": (276, 53.14474281170343) + "thoracic cavity": (40, 6.974341918899202), + "shell": (276, 52.878054197629744) } for name in expectedSizes3d: term = get_body_term(name) @@ -116,12 +116,12 @@ def test_wholebody2_core(self): expectedSizes2d = { "abdominal cavity boundary": (64, 27.428203203724674), "diaphragm": (20, 3.0778936559347208), - "left arm skin epidermis": (68, 22.627795339108015), + "left arm skin epidermis": (68, 22.433540462588258), "left leg skin epidermis": (68, 55.24949927903832), - "right arm skin epidermis": (68, 22.62779533911023), + "right arm skin epidermis": (68, 22.433540462588258), "right leg skin epidermis": (68, 55.24949927903832), - "skin epidermis": (388, 229.91482811362502), - "thoracic cavity boundary": (64, 20.606925296069125) + "skin epidermis": (388, 228.9017722286146), + "thoracic cavity boundary": (64, 20.2627556651649) } for name in expectedSizes2d: term = get_body_term(name) @@ -137,7 +137,7 @@ def test_wholebody2_core(self): self.assertAlmostEqual(surfaceArea, expectedSizes2d[name][1], delta=tol) expectedSizes1d = { - "spinal cord": (8, 10.854535713176467) + "spinal cord": (8, 10.85369421002332) } for name in expectedSizes1d: term = get_body_term(name) @@ -184,8 +184,8 @@ def test_wholebody2_tube(self): self.assertTrue(coordinates.isValid()) minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) tol = 1.0E-4 - assertAlmostEqualList(self, minimums, [0.0, -3.7000751482231564, -1.25], tol) - assertAlmostEqualList(self, maximums, [20.48332368641937, 3.7000751482231564, 2.15], tol) + assertAlmostEqualList(self, minimums, [0.0, -3.650833433150559, -1.25], tol) + assertAlmostEqualList(self, maximums, [20.48332368641937, 3.650833433150559, 2.15], tol) with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) @@ -212,13 +212,13 @@ def test_wholebody2_tube(self): result, innerSurfaceArea = innerSurfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 53.1439368556073, delta=tol) - self.assertAlmostEqual(outerSurfaceArea, 225.80985677023105, delta=tol) - self.assertAlmostEqual(innerSurfaceArea, 155.85891324174113, delta=tol) + self.assertAlmostEqual(volume, 52.87781598884186, delta=tol) + self.assertAlmostEqual(outerSurfaceArea, 224.9456647093771, delta=tol) + self.assertAlmostEqual(innerSurfaceArea, 155.4114878443358, delta=tol) # check some annotationGroups: expectedSizes2d = { - "skin epidermis": (320, 229.05761249752587) + "skin epidermis": (320, 228.11749988635236) } for name in expectedSizes2d: term = get_body_term(name)