ANCFThinPlateRevisionTest.py

You can view and download this file on Github: ANCFThinPlateRevisionTest.py

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Tests of ObjectANCFThinPlate and exudyn.shells (#2857, #2858): the integration modes against
#           F L^3/(3 EI), the material curvature measure, large-rotation pure bending, the first
#           eigenfrequency, the mass, gravity, variable thickness (4 and 12 values), Kelvin-Voigt
#           damping (no damping of a rigid motion, the rule of the bending damping), the output
#           variables, ShellMesh with a surface map and a thickness function, and the visualization
#           settings of shells. Small meshes and linear or few-step static solutions; damping is checked
#           with the linearized system. Each test prints the measured and the reference value.
#           The result is the number of tests that passed.
#
# Author:   Michael Pieber, Johannes Gerstmayr
# Date:     2026-10-01, ported 2026-10-05
#
# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

import io
import unittest
import numpy as np
import exudyn as exu
from exudyn.utilities import (MarkerNodePosition, LoadForceVector, MarkerNodeCoordinate, NodePointGround,
                              CoordinateConstraint)
from exudyn.shells import ShellMesh, SurfaceMap, AddNodeConstraints
from exudyn.utilities import MarkerBodyMass, LoadMassProportional

#material and strip geometry (Poisson ratio 0: plate strip = Euler-Bernoulli beam)
youngsModulus = 2.1e11
poissonRatio = 0.
density = 7800.
stripLength = 1.0
stripWidth = 0.1
stripThickness = 0.01
bendingStiffness = youngsModulus * stripWidth * stripThickness**3 / 12.     #EI of the strip
plateBendingStiffness = youngsModulus * stripThickness**3 / 12.               #D per unit width (nu = 0)
tipForce = 1.
gravity = 9.81
numberOfPlateDOFsPerNode = 9                   #position, slope x, slope y
firstBeamEigenvalueFactor = 1.87510407         #lambda_1 of the clamped-free Euler-Bernoulli beam
numberOfTipLoadSteps = 1
numberOfBendingLoadSteps = 5


def StaticSettings(numberOfLoadSteps=numberOfTipLoadSteps):
    """Return simulation settings for the static tests (no solution file, sparse solver).
    @param numberOfLoadSteps: number of static load steps
    @return: exu.SimulationSettings"""
    simulationSettings = exu.SimulationSettings()
    simulationSettings.staticSolver.numberOfLoadSteps = numberOfLoadSteps
    simulationSettings.staticSolver.newton.relativeTolerance = 1e-10
    simulationSettings.staticSolver.newton.absoluteTolerance = 1e-3   #round-off floor of the residual (N)
    simulationSettings.staticSolver.verboseMode = 0
    simulationSettings.solution.file.write = False
    simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
    return simulationSettings


def CreateStrip(numberOfElements=4, integrationMode=0, alongY=False, clamp=True, thicknessField=None,
                massProportionalLoad=None, stiffnessProportionalDamping=0.):
    """Create a plate strip (stripLength x stripWidth, one element across the width) with ShellMesh, optionally
    clamped at the edge at the origin.
    @param numberOfElements: number of elements along the strip
    @param integrationMode: useReducedOrderIntegration of the elements
    @param alongY: True: strip axis along global y, False: along global x
    @param clamp: True: all 9 coordinates of the nodes at the origin edge are fixed
    @param thicknessField: function f(x, y) of the thickness, evaluated at the nodes (4 values per element)
    @param massProportionalLoad: a load per mass, e.g. gravity, on every element
    @param stiffnessProportionalDamping: Kelvin-Voigt damping coefficient of ShellMesh
    @return: (system container, mbs, dictionary of nodes, elements, edgeNodeNumbers and nodeId, list of tip node numbers)"""
    systemContainer = exu.SystemContainer()
    mbs = systemContainer.AddSystem()
    (lengthX, lengthY, nx, ny) = ((stripWidth, stripLength, 1, numberOfElements) if alongY
                                  else (stripLength, stripWidth, numberOfElements, 1))
    plate = ShellMesh(vertices=[[0, 0, 0], [lengthX, 0, 0], [lengthX, lengthY, 0], [0, lengthY, 0]],
                      numberOfElementsX=nx, numberOfElementsY=ny, youngsModulus=youngsModulus,
                      poissonsRatio=poissonRatio, density=density, thickness=stripThickness,
                      stiffnessProportionalDamping=stiffnessProportionalDamping)
    if thicknessField is not None:
        plate.thicknessAtNodes = [thicknessField(ix * lengthX / nx, iy * lengthY / ny)
                                  for iy in range(ny + 1) for ix in range(nx + 1)]
    plate.CreateANCFThinPlateElements(mbs, useReducedOrderIntegration=integrationMode)
    if massProportionalLoad is not None:
        for elementNumber in plate.elementNumbers:
            mbs.AddLoad(LoadMassProportional(markerNumber=mbs.AddMarker(MarkerBodyMass(bodyNumber=elementNumber)),
                                             loadVector=massProportionalLoad))
    built = {'nodes': plate.nodeNumbers, 'elements': plate.elementNumbers,
             'edgeNodeNumbers': plate.boundaryNodeNumbers, 'nodeId': lambda iy, ix: iy * (nx + 1) + ix}
    if alongY:
        clampNodes, tipNodes = built['edgeNodeNumbers']['bottom'], built['edgeNodeNumbers']['top']
    else:
        clampNodes, tipNodes = built['edgeNodeNumbers']['left'], built['edgeNodeNumbers']['right']
    if clamp:
        for nodeNumber in clampNodes:
            AddNodeConstraints(mbs, nodeNumber, dofs=list(range(numberOfPlateDOFsPerNode)))
    return systemContainer, mbs, built, tipNodes


def AddTipLoad(mbs, tipNodes, force=tipForce):
    """Distribute a force in -z equally on the tip nodes.
    @param mbs: main system
    @param tipNodes: list of tip node numbers
    @param force: total force (N)"""
    for nodeNumber in tipNodes:
        markerNumber = mbs.AddMarker(MarkerNodePosition(nodeNumber=nodeNumber))
        mbs.AddLoad(LoadForceVector(markerNumber=markerNumber, loadVector=[0, 0, -force / len(tipNodes)]))


def TipDeflection(mbs, tipNodes):
    """Return the mean downward (-z) displacement of the tip nodes.
    @param mbs: main system (solved)
    @param tipNodes: list of tip node numbers
    @return: deflection (m), positive downward"""
    return -np.mean([mbs.GetNodeOutput(nodeNumber, exu.OutputVariableType.Displacement)[2] for nodeNumber in tipNodes])


def RigidMotionVectors(mbs, rotationMatrix, angularVelocity):
    """Coordinates of a rigid rotation of all plate nodes and the velocities of a rigid angular velocity in
    this rotated configuration (positions and both slope vectors rotate with the body).
    @param mbs: main system (assembled, NodePointSlope12 nodes only)
    @param rotationMatrix: 3x3 rotation of the reference configuration
    @param angularVelocity: angular velocity vector (rad/s)
    @return: (ODE2 coordinates, ODE2 velocities) as numpy arrays"""
    coordinates = np.zeros(mbs.systemData.ODE2Size())
    velocities = np.zeros(mbs.systemData.ODE2Size())
    for nodeNumber in range(mbs.systemData.NumberOfNodes()):
        referenceCoordinates = np.array(mbs.GetNodeParameter(nodeNumber, 'referenceCoordinates'))
        firstIndex = mbs.GetNodeODE2Index(nodeNumber)
        for vectorIndex in range(3):
            referenceVector = referenceCoordinates[3 * vectorIndex:3 * vectorIndex + 3]
            rotatedVector = rotationMatrix @ referenceVector
            coordinates[firstIndex + 3 * vectorIndex:firstIndex + 3 * vectorIndex + 3] = rotatedVector - referenceVector
            velocities[firstIndex + 3 * vectorIndex:firstIndex + 3 * vectorIndex + 3] = np.cross(angularVelocity, rotatedVector)
    return coordinates, velocities


def PrintComparison(label, measured, reference):
    """Verbose test output: measured value, reference value and relative deviation.
    @param label: description of the quantity
    @param measured: measured value
    @param reference: reference value (nonzero)"""
    print(f'\n    {label}: measured {measured:.9e}, reference {reference:.9e}, '
          f'relative deviation {measured / reference - 1:+.3e}', end='')


class TestANCFThinPlateElement(unittest.TestCase):
    """Kinematics, integration, mass, gravity and thickness of ObjectANCFThinPlate."""

    def testCantileverTipDeflectionAllIntegrationModes(self):
        """Linear cantilever with tip load: cubic Hermite elements are exact, F L^3/(3 EI), for modes 0, 1, 2;
        mode 2 is an alias of mode 1 (bitwise identical)."""
        deflections = {}
        for integrationMode in [0, 1, 2]:
            systemContainer, mbs, built, tipNodes = CreateStrip(numberOfElements=2, integrationMode=integrationMode)
            AddTipLoad(mbs, tipNodes)
            mbs.Assemble()
            mbs.SolveStatic(StaticSettings())
            deflections[integrationMode] = TipDeflection(mbs, tipNodes)
            reference = tipForce * stripLength**3 / (3 * bendingStiffness)
            PrintComparison(f'tip deflection mode {integrationMode}', deflections[integrationMode], reference)
            self.assertAlmostEqual(deflections[integrationMode] / reference, 1., delta=1e-6)
        self.assertEqual(deflections[1], deflections[2])

    @unittest.skipIf('[FAST]' in exu.config.Version(True),
                     'exudynCPPfast compiles the error checks away by design (#2911)')
    def testInvalidIntegrationModeRaises(self):
        """useReducedOrderIntegration = 3 must stop the solver with an error."""
        systemContainer, mbs, built, tipNodes = CreateStrip(numberOfElements=2, integrationMode=3)
        mbs.Assemble()
        printToConsole = exu.config.printToConsole
        exu.config.printToConsole = False    #the expected error message would clutter the test output
        try:
            with self.assertRaises(ValueError):
                mbs.SolveStatic(StaticSettings())
        finally:
            exu.config.printToConsole = printToConsole

    def testMaterialCurvatureOfStretchedArc(self):
        """Material curvature measure: a strip placed exactly on a circular arc of radius rho and stretched by
        lambda has |kappa_xx| = lambda/rho (geometric measure: 1/rho) and Green-Lagrange strain (lambda^2-1)/2.
        Pure kinematics, no solver."""
        stretch, radius = 1.2, 1.0
        systemContainer, mbs, built, tipNodes = CreateStrip(numberOfElements=4, clamp=False)
        mbs.Assemble()
        coordinates = np.zeros(mbs.systemData.ODE2Size())
        for nodeNumber in range(mbs.systemData.NumberOfNodes()):
            referenceCoordinates = np.array(mbs.GetNodeParameter(nodeNumber, 'referenceCoordinates'))
            arcAngle = stretch * referenceCoordinates[0] / radius
            position = np.array([radius * np.sin(arcAngle), referenceCoordinates[1], radius * (1 - np.cos(arcAngle))])
            slopeX = stretch * np.array([np.cos(arcAngle), 0, np.sin(arcAngle)])
            firstIndex = mbs.GetNodeODE2Index(nodeNumber)
            coordinates[firstIndex:firstIndex + 3] = position - referenceCoordinates[0:3]
            coordinates[firstIndex + 3:firstIndex + 6] = slopeX - referenceCoordinates[3:6]
        mbs.systemData.SetODE2Coordinates(coordinates.tolist(), configuration=exu.ConfigurationType.Current)
        for elementNumber in built['elements']:
            curvature = mbs.GetObjectOutputBody(elementNumber, exu.OutputVariableType.CurvatureLocal,
                                                localPosition=[0, 0, 0], configuration=exu.ConfigurationType.Current)
            strain = mbs.GetObjectOutputBody(elementNumber, exu.OutputVariableType.StrainLocal,
                                             localPosition=[0, 0, 0], configuration=exu.ConfigurationType.Current)
            PrintComparison(f'element {elementNumber} |kappa_xx|', abs(curvature[0]), stretch / radius)
            self.assertAlmostEqual(abs(curvature[0]) / (stretch / radius), 1., delta=1e-2)   #interpolation error 0.4 %
            self.assertAlmostEqual(strain[0], 0.5 * (stretch**2 - 1), delta=1e-4)
            self.assertAlmostEqual(curvature[1], 0., delta=1e-10)

    def testPureBendingLargeRotation(self):
        """Tip slope rotated by 1 rad: the strip forms a circular arc of radius L/theta; curvature theta/L,
        membrane strain small compared with the surface bending strain."""
        rotationAngle = 1.0
        systemContainer, mbs, built, tipNodes = CreateStrip(numberOfElements=4)
        groundNode = mbs.AddNode(NodePointGround())
        targetSlope = [np.cos(rotationAngle), 0., np.sin(rotationAngle)]
        for nodeNumber in tipNodes:
            referenceCoordinates = mbs.GetNodeParameter(nodeNumber, 'referenceCoordinates')
            for component in range(3):
                groundMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=groundNode, coordinate=0))
                slopeMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodeNumber, coordinate=3 + component))
                mbs.AddObject(CoordinateConstraint(markerNumbers=[groundMarker, slopeMarker],
                                                   offset=targetSlope[component] - referenceCoordinates[3 + component]))
        mbs.Assemble()
        mbs.SolveStatic(StaticSettings(numberOfBendingLoadSteps))
        tipPosition = np.mean([mbs.GetNodeOutput(nodeNumber, exu.OutputVariableType.Position) for nodeNumber in tipNodes], axis=0)
        exactTipPosition = stripLength / rotationAngle * np.array([np.sin(rotationAngle), 1 - np.cos(rotationAngle)])
        arcError = np.linalg.norm(tipPosition[[0, 2]] - exactTipPosition)
        print(f'\n    tip position error {arcError:.2e} m', end='')
        self.assertLess(arcError, 1e-5)
        for elementNumber in built['elements']:
            curvature = mbs.GetObjectOutputBody(elementNumber, exu.OutputVariableType.CurvatureLocal, localPosition=[0, 0, 0])
            membraneStrain = mbs.GetObjectOutputBody(elementNumber, exu.OutputVariableType.StrainLocal, localPosition=[0, 0, 0])
            surfaceStrain = mbs.GetObjectOutputBody(elementNumber, exu.OutputVariableType.StrainLocal, localPosition=[0, 0, 1])
            PrintComparison(f'element {elementNumber} |kappa_xx|', abs(curvature[0]), rotationAngle / stripLength)
            self.assertAlmostEqual(abs(curvature[0]) / (rotationAngle / stripLength), 1., delta=5e-3)
            self.assertLess(abs(membraneStrain[0] / surfaceStrain[0]), 2e-3)

    def testFirstEigenfrequencyCantilever(self):
        """First bending eigenfrequency of the clamped strip vs. Euler-Bernoulli (mass and stiffness)."""
        systemContainer, mbs, built, tipNodes = CreateStrip(numberOfElements=4, clamp=False)
        mbs.Assemble()
        constrainedCoordinates = []
        for nodeNumber in built['edgeNodeNumbers']['left']:
            firstIndex = mbs.GetNodeODE2Index(nodeNumber)
            constrainedCoordinates += list(range(firstIndex, firstIndex + numberOfPlateDOFsPerNode))
        frequencies = mbs.ComputeODE2Eigenvalues(constrainedCoordinates=constrainedCoordinates,
                                                 convert2Frequencies=True, useSparseSolver=False)[0]
        reference = (firstBeamEigenvalueFactor**2 / (2 * np.pi)
                     * np.sqrt(bendingStiffness / (density * stripWidth * stripThickness * stripLength**4)))
        PrintComparison('first eigenfrequency (Hz)', frequencies[0], reference)
        self.assertAlmostEqual(frequencies[0] / reference, 1., delta=1e-4)

    def testMassOfLinearlyTaperedPlate(self):
        """Mass matrix with thickness h(x) = h0 (1 - x/(2L)): total mass rho b h0 L 3/4 (exact)."""
        systemContainer, mbs, built, tipNodes = CreateStrip(
            numberOfElements=4, clamp=False, thicknessField=lambda x, y: stripThickness * (1 - 0.5 * x / stripLength))
        mbs.Assemble()
        massMatrix = mbs.ComputeLinearizedSystem(StaticSettings())[0]
        translationZ = np.zeros(mbs.systemData.ODE2Size())
        for nodeNumber in range(mbs.systemData.NumberOfNodes()):
            translationZ[mbs.GetNodeODE2Index(nodeNumber) + 2] = 1.
        reference = density * stripWidth * stripThickness * stripLength * 0.75
        PrintComparison('total mass (kg)', translationZ @ massMatrix @ translationZ, reference)
        self.assertAlmostEqual(translationZ @ massMatrix @ translationZ / reference, 1., delta=1e-10)

    def testGravityTipDeflection(self):
        """Self-weight (LoadMassProportional on each element): tip deflection q L^4/(8 EI)."""
        systemContainer, mbs, built, tipNodes = CreateStrip(numberOfElements=4, massProportionalLoad=[0, 0, -gravity])
        mbs.Assemble()
        mbs.SolveStatic(StaticSettings())
        lineLoad = density * stripWidth * stripThickness * gravity
        reference = lineLoad * stripLength**4 / (8 * bendingStiffness)
        PrintComparison('tip deflection under self-weight', TipDeflection(mbs, tipNodes), reference)
        self.assertAlmostEqual(TipDeflection(mbs, tipNodes) / reference, 1., delta=1e-4)

    def testUniformVariableThicknessEqualsConstant(self):
        """4 equal nodal thicknesses (point-wise stiffness D_base h) give the same result as constant thickness."""
        deflections = []
        for builderOptions in [{}, {'thicknessField': lambda x, y: stripThickness}]:
            systemContainer, mbs, built, tipNodes = CreateStrip(numberOfElements=2, **builderOptions)
            AddTipLoad(mbs, tipNodes)
            mbs.Assemble()
            mbs.SolveStatic(StaticSettings())
            deflections.append(TipDeflection(mbs, tipNodes))
        PrintComparison('tip deflection, 4 equal thickness values', deflections[1], deflections[0])
        self.assertAlmostEqual(deflections[1] / deflections[0], 1., delta=1e-10)

    def testLinearTaperTipDeflection(self):
        """Thickness h(x) = h0 (1 - x/(2L)), bending stiffness point-wise with h^3: tip deflection
        integral of F (L-x)^2 / (E b h(x)^3/12) dx."""
        thicknessFunction = lambda x, y: stripThickness * (1 - 0.5 * x / stripLength)
        systemContainer, mbs, built, tipNodes = CreateStrip(numberOfElements=4, thicknessField=thicknessFunction)
        AddTipLoad(mbs, tipNodes)
        mbs.Assemble()
        mbs.SolveStatic(StaticSettings())
        positions = np.linspace(0, stripLength, 20001)
        integrand = tipForce * (stripLength - positions)**2 / (youngsModulus * stripWidth * thicknessFunction(positions, 0)**3 / 12)
        reference = np.sum(0.5 * (integrand[1:] + integrand[:-1]) * np.diff(positions))   #trapezoidal rule
        PrintComparison('tip deflection, linear taper', TipDeflection(mbs, tipNodes), reference)
        self.assertAlmostEqual(TipDeflection(mbs, tipNodes) / reference, 1., delta=1e-3)

    def testTwelveThicknessValuesEqualConstant(self):
        """12 thickness values [h_i, h_s,i, h_t,i] with zero gradients give the deflection of the constant thickness."""
        deflections = []
        for twelveValues in [False, True]:
            systemContainer, mbs, built, tipNodes = CreateStrip(numberOfElements=2)
            if twelveValues:
                for elementNumber in built['elements']:
                    mbs.SetObjectParameter(elementNumber, 'thickness', [stripThickness, 0., 0.] * 4)
            AddTipLoad(mbs, tipNodes)
            mbs.Assemble()
            mbs.SolveStatic(StaticSettings())
            deflections.append(TipDeflection(mbs, tipNodes))
        PrintComparison('tip deflection, 12 thickness values', deflections[1], deflections[0])
        self.assertAlmostEqual(deflections[1] / deflections[0], 1., delta=1e-12)

class TestANCFThinPlateDamping(unittest.TestCase):
    """Kelvin-Voigt damping, checked with the damping matrix of the linearized system."""

    def DampingMatrix(self, membraneDamping, bendingDamping, coordinates=None):
        """Damping matrix of a free 2-element strip.
        @param membraneDamping: stiffnessProportionalDamping (s)
        @param bendingDamping: bendingStiffnessProportionalDamping (s), negative: use membraneDamping
        @param coordinates: optional ODE2 coordinates (initial and current configuration)
        @return: (mbs, stiffness matrix, damping matrix)"""
        systemContainer, mbs, built, tipNodes = CreateStrip(numberOfElements=2, clamp=False,
                                                            stiffnessProportionalDamping=membraneDamping)
        for elementNumber in built['elements']:
            mbs.SetObjectParameter(elementNumber, 'bendingStiffnessProportionalDamping', bendingDamping)
        mbs.Assemble()
        if coordinates is not None:
            for configuration in [exu.ConfigurationType.Initial, exu.ConfigurationType.Current]:
                mbs.systemData.SetODE2Coordinates(coordinates.tolist(), configuration=configuration)
        massMatrix, stiffnessMatrix, dampingMatrix = mbs.ComputeLinearizedSystem(StaticSettings())
        return mbs, stiffnessMatrix, dampingMatrix

    def testNoDampingOfRigidRotation(self):
        """In a finitely rotated configuration, a rigid angular velocity produces no damping force (strain
        rates vanish), while a bending velocity does; also no elastic force change for rigid rotation."""
        rotationAngle = 1.0
        rotationMatrix = np.array([[np.cos(rotationAngle), 0, np.sin(rotationAngle)], [0, 1, 0],
                                   [-np.sin(rotationAngle), 0, np.cos(rotationAngle)]])
        mbs, stiffnessMatrix, dampingMatrix = self.DampingMatrix(1e-3, -1.)
        coordinates, rigidVelocities = RigidMotionVectors(mbs, rotationMatrix, np.array([0.3, -0.5, 0.7]))
        mbs, stiffnessMatrix, dampingMatrix = self.DampingMatrix(1e-3, -1., coordinates)
        bendingVelocities = np.zeros(mbs.systemData.ODE2Size())
        for nodeNumber in range(mbs.systemData.NumberOfNodes()):
            referenceCoordinates = np.array(mbs.GetNodeParameter(nodeNumber, 'referenceCoordinates'))
            firstIndex = mbs.GetNodeODE2Index(nodeNumber)
            bendingVelocities[firstIndex:firstIndex + 3] = rotationMatrix @ np.array([0, 0, referenceCoordinates[0]**2])
        rigidForce = np.linalg.norm(dampingMatrix @ rigidVelocities)
        bendingForce = np.linalg.norm(dampingMatrix @ bendingVelocities)
        print(f'\n    |D v_rigid| = {rigidForce:.3e}, |D v_bending| = {bendingForce:.3e}, '
              f'|K v_rigid|/|K v_bending| = {np.linalg.norm(stiffnessMatrix @ rigidVelocities) / np.linalg.norm(stiffnessMatrix @ bendingVelocities):.3e}', end='')
        self.assertGreater(bendingForce, 1.)
        self.assertLess(rigidForce / bendingForce, 1e-9)
        self.assertLess(np.linalg.norm(stiffnessMatrix @ rigidVelocities) / np.linalg.norm(stiffnessMatrix @ bendingVelocities), 1e-9)

    def testBendingDampingParameterRule(self):
        """Flat strip: out-of-plane velocity w = x^2 only strains in bending, in-plane velocity u = x only in the
        membrane. beta_kappa >= 0 is used as is (0 switches bending damping off), negative uses beta_epsilon."""
        cases = [(1e-3, 0., False, True), (1e-3, -1., True, True), (0., 1e-3, True, False), (0., -1., False, False)]
        for membraneDamping, bendingDamping, bendingDamped, membraneDamped in cases:
            mbs, stiffnessMatrix, dampingMatrix = self.DampingMatrix(membraneDamping, bendingDamping)
            bendingVelocities = np.zeros(mbs.systemData.ODE2Size())
            stretchVelocities = np.zeros(mbs.systemData.ODE2Size())
            for nodeNumber in range(mbs.systemData.NumberOfNodes()):
                referenceCoordinates = np.array(mbs.GetNodeParameter(nodeNumber, 'referenceCoordinates'))
                firstIndex = mbs.GetNodeODE2Index(nodeNumber)
                bendingVelocities[firstIndex + 2] = referenceCoordinates[0]**2          #w = x^2
                bendingVelocities[firstIndex + 5] = 2 * referenceCoordinates[0]         #w_x = 2x
                stretchVelocities[firstIndex] = referenceCoordinates[0]                 #u = x
                stretchVelocities[firstIndex + 3] = 1.                                  #u_x = 1
            bendingForce = np.linalg.norm(dampingMatrix @ bendingVelocities)
            stretchForce = np.linalg.norm(dampingMatrix @ stretchVelocities)
            print(f'\n    beta_epsilon = {membraneDamping}, beta_kappa = {bendingDamping}: '
                  f'|D v_bending| = {bendingForce:.3e}, |D v_stretch| = {stretchForce:.3e}', end='')
            self.assertEqual(bendingForce > 1e-6, bendingDamped)
            self.assertEqual(stretchForce > 1e-6, membraneDamped)


class TestANCFThinPlateOutputs(unittest.TestCase):
    """Output variables StressLocal, StrainLocal, CurvatureLocal, ForceLocal, TorqueLocal."""

    def SolvedTipLoadStrip(self, alongY):
        """Clamped 4-element strip with tip load in -z (concave down), solved.
        @param alongY: strip axis along global y
        @return: (mbs, built dictionary)"""
        systemContainer, mbs, built, tipNodes = CreateStrip(numberOfElements=4, alongY=alongY)
        AddTipLoad(mbs, tipNodes)
        mbs.Assemble()
        mbs.SolveStatic(StaticSettings())
        return mbs, built

    def testOutputVectorSizes(self):
        """StressLocal and StrainLocal have 6 components, CurvatureLocal, ForceLocal and TorqueLocal 3."""
        mbs, built = self.SolvedTipLoadStrip(False)
        expectedSizes = {exu.OutputVariableType.StressLocal: 6, exu.OutputVariableType.StrainLocal: 6,
                         exu.OutputVariableType.CurvatureLocal: 3, exu.OutputVariableType.ForceLocal: 3,
                         exu.OutputVariableType.TorqueLocal: 3}
        for outputVariableType, expectedSize in expectedSizes.items():
            value = mbs.GetObjectOutputBody(built['elements'][0], outputVariableType, localPosition=[0, 0, 0])
            print(f'\n    {outputVariableType}: {len(value)} components', end='')
            self.assertEqual(len(value), expectedSize)

    def testStressStrainAndBendingMoment(self):
        """Top fiber (zeta = +1) near the clamp in tension for concave-down bending along x and along y;
        sigma = E epsilon (nu = 0); |sigma| = 6 m/h^2 with m = F (L - s)/b at the element center."""
        for alongY in [False, True]:
            mbs, built = self.SolvedTipLoadStrip(alongY)
            component = 1 if alongY else 0
            elementNumber = built['elements'][0]
            stressTop = mbs.GetObjectOutputBody(elementNumber, exu.OutputVariableType.StressLocal, localPosition=[0, 0, 1])[component]
            strainTop = mbs.GetObjectOutputBody(elementNumber, exu.OutputVariableType.StrainLocal, localPosition=[0, 0, 1])[component]
            momentPerWidth = tipForce * (stripLength - 0.5 * stripLength / 4) / stripWidth
            PrintComparison(f'strip along {"y" if alongY else "x"}: top stress', stressTop, 6 * momentPerWidth / stripThickness**2)
            self.assertGreater(stressTop, 0.)
            self.assertAlmostEqual(stressTop / (youngsModulus * strainTop), 1., delta=1e-8)
            self.assertAlmostEqual(stressTop / (6 * momentPerWidth / stripThickness**2), 1., delta=1e-5)   #nonlinear static solution

    def testMomentEqualsPlateStiffnessTimesCurvature(self):
        """TorqueLocal = E h^3/12 CurvatureLocal component-wise (nu = 0), i.e. moment and curvature use the same
        sign convention, for strips along x and along y."""
        for alongY in [False, True]:
            mbs, built = self.SolvedTipLoadStrip(alongY)
            component = 1 if alongY else 0
            curvature = mbs.GetObjectOutputBody(built['elements'][0], exu.OutputVariableType.CurvatureLocal, localPosition=[0, 0, 0])
            moment = mbs.GetObjectOutputBody(built['elements'][0], exu.OutputVariableType.TorqueLocal, localPosition=[0, 0, 0])
            PrintComparison(f'strip along {"y" if alongY else "x"}: moment/curvature', moment[component] / curvature[component],
                            plateBendingStiffness)
            self.assertAlmostEqual(moment[component] / curvature[component] / plateBendingStiffness, 1., delta=1e-8)

    def testCurvatureSignSameInBothDirections(self):
        """For the same concave-down bending, CurvatureLocal component 0 (strip along x) and component 1 (strip
        along y) have the same sign, the sign of ObjectANCFCable2D: negative (#2864)."""
        curvatures = []
        for alongY in [False, True]:
            mbs, built = self.SolvedTipLoadStrip(alongY)
            curvature = mbs.GetObjectOutputBody(built['elements'][0], exu.OutputVariableType.CurvatureLocal, localPosition=[0, 0, 0])
            curvatures.append(curvature[1 if alongY else 0])
        print(f'\n    curvature along x: {curvatures[0]:+.4e}, along y: {curvatures[1]:+.4e}', end='')
        self.assertEqual(np.sign(curvatures[0]), np.sign(curvatures[1]))
        self.assertLess(curvatures[0], 0.)


class TestShellsModule(unittest.TestCase):
    """exudyn.shells: ShellMesh, surface maps and thickness functions."""

    def testShellMeshTopology(self):
        """3 x 2 elements: 12 nodes, 6 elements, edge node lists, node positions on the grid."""
        systemContainer = exu.SystemContainer()
        mbs = systemContainer.AddSystem()
        numberOfElementsX, numberOfElementsY, lengthX, lengthY = 3, 2, 0.6, 0.4
        origin = np.array([0.1, -0.2, 0.3])
        plate = ShellMesh(vertices=[origin, origin + [lengthX, 0, 0], origin + [lengthX, lengthY, 0], origin + [0, lengthY, 0]],
                          numberOfElementsX=numberOfElementsX, numberOfElementsY=numberOfElementsY,
                          youngsModulus=youngsModulus, poissonsRatio=0.3, density=density, thickness=stripThickness)
        plate.CreateANCFThinPlateElements(mbs)
        print(f'\n    nodes {len(plate.nodeNumbers)}, elements {len(plate.elementNumbers)}, edge nodes '
              f'{ {key: len(value) for key, value in plate.boundaryNodeNumbers.items()} }', end='')
        self.assertEqual(len(plate.nodeNumbers), (numberOfElementsX + 1) * (numberOfElementsY + 1))
        self.assertEqual(len(plate.elementNumbers), numberOfElementsX * numberOfElementsY)
        self.assertEqual(len(plate.boundaryNodeNumbers['left']), numberOfElementsY + 1)
        self.assertEqual(len(plate.boundaryNodeNumbers['bottom']), numberOfElementsX + 1)
        for indexY in range(numberOfElementsY + 1):
            for indexX in range(numberOfElementsX + 1):
                nodeNumber = plate.nodeNumbers[indexY * (numberOfElementsX + 1) + indexX]
                position = np.array(mbs.GetNodeParameter(nodeNumber, 'referenceCoordinates'))[0:3]
                expected = origin + [indexX * lengthX / numberOfElementsX, indexY * lengthY / numberOfElementsY, 0]
                self.assertLess(np.linalg.norm(position - expected), 1e-12)

    def testHemisphericalMapNodesOnSphere(self):
        """SurfaceMap.HemisphericalShell: all nodes on the sphere of radius R, both slope vectors tangent to it."""
        systemContainer = exu.SystemContainer()
        mbs = systemContainer.AddSystem()
        radius = 2.0
        plate = ShellMesh(vertices=[[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]], numberOfElementsX=4, numberOfElementsY=4,
                          youngsModulus=youngsModulus, poissonsRatio=0.3, density=density, thickness=stripThickness,
                          surfaceMap=SurfaceMap.HemisphericalShell(1., 1., radius, np.pi / 3))
        plate.CreateANCFThinPlateElements(mbs)
        radiusError, normalSlopeComponent = 0., 0.
        for nodeNumber in plate.nodeNumbers:
            referenceCoordinates = np.array(mbs.GetNodeParameter(nodeNumber, 'referenceCoordinates'))
            position = referenceCoordinates[0:3]
            radiusError = max(radiusError, abs(np.linalg.norm(position) - radius))
            for slope in [referenceCoordinates[3:6], referenceCoordinates[6:9]]:
                if np.linalg.norm(slope) > 0:       #degenerate slope at the pole
                    normalSlopeComponent = max(normalSlopeComponent, abs(slope @ position) / (np.linalg.norm(slope) * radius))
        print(f'\n    radius error {radiusError:.2e}, normal slope component {normalSlopeComponent:.2e}', end='')
        self.assertLess(radiusError, 1e-12)
        self.assertLess(normalSlopeComponent, 1e-12)

    def testThicknessFunctionOfARotatedStrip(self):
        """A strip turned by 30 degrees in the x-y plane with the thickness h0 (1 - s/(2L)) along its axis s, given as a
        function of the global x and y: the gradients along the node slopes give 12 thickness values per element, and
        the tip deflection is the one of the tapered cantilever (#2864)."""
        angle = np.pi / 6
        axis, across = np.array([np.cos(angle), np.sin(angle), 0.]), np.array([-np.sin(angle), np.cos(angle), 0.])
        def ThicknessFunction(x, y):
            return stripThickness * (1 - 0.5 * (x * axis[0] + y * axis[1]) / stripLength)
        systemContainer = exu.SystemContainer()
        mbs = systemContainer.AddSystem()
        plate = ShellMesh(vertices=[[0, 0, 0], stripLength * axis, stripLength * axis + stripWidth * across, stripWidth * across],
                          numberOfElementsX=8, numberOfElementsY=1, youngsModulus=youngsModulus, poissonsRatio=poissonRatio,
                          density=density, thickness=stripThickness, thicknessFunction=ThicknessFunction)
        plate.CreateANCFThinPlateElements(mbs)
        self.assertEqual(len(mbs.GetObjectParameter(plate.elementNumbers[0], 'thickness')), 12)
        for nodeNumber in plate.boundaryNodeNumbers['left']:
            AddNodeConstraints(mbs, nodeNumber, dofs=list(range(numberOfPlateDOFsPerNode)))
        tipNodes = plate.boundaryNodeNumbers['right']
        for nodeNumber in tipNodes:
            mbs.AddLoad(LoadForceVector(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=nodeNumber)),
                                        loadVector=[0, 0, -tipForce / len(tipNodes)]))
        mbs.Assemble()
        mbs.SolveStatic(StaticSettings())
        positions = np.linspace(0, stripLength, 20001)
        thickness = stripThickness * (1 - 0.5 * positions / stripLength)
        integrand = tipForce * (stripLength - positions)**2 / (youngsModulus * stripWidth * thickness**3 / 12)
        reference = np.sum(0.5 * (integrand[1:] + integrand[:-1]) * np.diff(positions))   #trapezoidal rule
        PrintComparison('tip deflection, thickness function of a rotated strip', TipDeflection(mbs, tipNodes), reference)
        self.assertAlmostEqual(TipDeflection(mbs, tipNodes) / reference, 1., delta=1e-3)


class TestVisualizationSettings(unittest.TestCase):
    """The settings of visualizationSettings.bodies.shells and their defaults."""

    def testShellVisualizationDefaults(self):
        """contourZeta, reducedInterpolation, drawNormal, drawNormalFactor, drawNormalColor, drawNormalLines."""
        shellSettings = exu.SystemContainer().visualizationSettings.bodies.shells
        self.assertEqual(shellSettings.contourZeta, 1.)
        self.assertTrue(shellSettings.reducedInterpolation)
        self.assertFalse(shellSettings.drawNormal)
        self.assertEqual(shellSettings.drawNormalFactor, 1.)
        self.assertTrue(np.allclose(list(shellSettings.drawNormalColor), [0.2, 0.2, 0.2, 1.]))
        self.assertTrue(shellSettings.drawNormalLines)


#run all tests; the result is the number of tests that passed
testIsActive = exu.sys.get('testIsActive', False)
result = unittest.TextTestRunner(verbosity=2 if not testIsActive else 0,
                                 stream=io.StringIO() if testIsActive else None).run(
    unittest.TestSuite(unittest.defaultTestLoader.loadTestsFromTestCase(testCase) for testCase in
                       [TestANCFThinPlateElement, TestANCFThinPlateDamping, TestANCFThinPlateOutputs, TestShellsModule,
                        TestVisualizationSettings]))
for (test, trace) in result.failures + result.errors:
    exu.Print('FAILED: ' + str(test) + '\n' + trace)
testResult = result.testsRun - len(result.failures) - len(result.errors) - len(result.unexpectedSuccesses)
exu.Print('solution of ANCFThinPlateRevisionTest=', testResult)
exu.sys['testResult'] = testResult