Coordinate System Transformation and implementation

Question

I am attempting to read analysis results from an output file generated by OOFEM (analysis engine) and then to apply the deflections read from the file onto my model using a combination of C# and VB.Net. I'm using devDept.Eyeshot for visualisation.

I am running into an issue when attempting to convert the OOFEM beam3d element output into the Global Coordinate system in my app. The OOFEM documentation on the matter:

Beam element for 3D linear analysis, based on Timchshenko hypothesis. The internal condensation of arbitrary DOF is supported and is performed in local system. On output, the local end-displacement and local end-forces are printed in local coordinate system. Requires the local system to be chosen according to main central axes of inertia. Local element coordinate system is determined by the following rules: l. let first element node has following coordinates (Xi, Yi, Zi,) and the second one (Xj, Yj, Zj), 2. direction vector of local X-axis is then a1 = (Xj – Xi, Yj – Yi, Zj – Zi), 3. local Y-axis direction vector lies in plane by local X-axis direction vector (a1) and given point (k-node with coordinates (Xk, Yk, Zk)) - so called reference node, 4. local Z-axis is then determined as vector product of local X-axis direction vector (a1) by vector (Xk – Xi, Yk – Yi, Zk – Zi), 5. local Y-axis is then determined as vector product of local Z-axis direction vector by local X-axis direction vector.

My code consists of a C# class for the CrossSection objects, which attempts to determine the KNode and the beam's X direction using the Cross Section properties. The Cross Sections are always defined as 2D region on the XY plane. Herewith an extract on the section calcs:

AreaProperties areaProperties = new AreaProperties();
List<Point3D> lstp3d = new List<Point3D>();
foreach (Point2D pt in coords){
    lstp3d.Add(new Point3D(pt.X, pt.Y, 0));
}
areaProperties.Add(lstp3d);
Centroid = new Point3D(Math.Round(areaProperties.Centroid.X, 5), Math.Round(areaProperties.Centroid.Y, 5), Math.Round(areaProperties.Centroid.Z, 5));
KNode = new Point3D(Centroid.X + width / 2, Centroid.Y, Centroid.Z);

Then in my own Beam3D class, I have:

public void SetLocalVectors(Node nde1, Node nde2)
{
    Vector3D localX = new Vector3D(nde1.ToPoint3D(), nde2.ToPoint3D());
    localX.Normalize();

    // Use the KNode property to get the direction to the k-node
    Vector3D toKNode = new Vector3D(nde1.ToPoint3D(), CrossSection.KNode);

    // Compute the Local Z-Axis
    Vector3D localZ = Vector3D.Cross(toKNode, localX);
    localZ.Normalize();

    // Compute the Local Y-Axis
    Vector3D localY = Vector3D.Cross(localZ, localX);
    localY.Normalize();

    LocalXDirection = localX;
    LocalYDirection = localY;
    LocalZDirection = localZ;
}

Where I attempt to set the beam’s local coordinate system so I can interpret the results correctly later, and apply to the global CS.

I’m reading the OOFEM results as follow:

public void ReadBeamOutput(string filePath, OOFEMProject project)
{
    if (File.Exists(filePath))
    {
        using (StreamReader reader = new StreamReader(filePath))
        {
            // Skip the header line
            reader.ReadLine();

            string line;
            while ((line = reader.ReadLine()) != null)
            {
                if (line.Contains("endStep"))
                {
                    break;
                }
                string[] parts = line.Split(';');

                double time = double.Parse(parts[0], CultureInfo.InvariantCulture);
                int beamId = int.Parse(parts[1]);
                double distanceFromI = Math.Round(double.Parse(parts[2], CultureInfo.InvariantCulture), 5);

                Beam3D beam = FindBeamById(beamId);
                if (beam != null)
                {
                    // Create a transformation matrix using the local axes
                    Matrix3D transformationMatrix = new Matrix3D(
                        beam.LocalXDirection.X, beam.LocalXDirection.Y, beam.LocalXDirection.Z, 0,
                        beam.LocalYDirection.X, beam.LocalYDirection.Y, beam.LocalYDirection.Z, 0,
                        beam.LocalZDirection.X, beam.LocalZDirection.Y, beam.LocalZDirection.Z, 0,
                        0, 0, 0, 1
                    );

                    // Convert local results to global coordinates
                    System.Windows.Media.Media3D.Vector3D localDeflection = new System.Windows.Media.Media3D.Vector3D(
                        Math.Round(double.Parse(parts[9], CultureInfo.InvariantCulture), 10),
                        Math.Round(double.Parse(parts[10], CultureInfo.InvariantCulture), 10),
                        Math.Round(double.Parse(parts[11], CultureInfo.InvariantCulture), 10)
                    );
                    System.Windows.Media.Media3D.Vector3D globalDeflection = System.Windows.Media.Media3D.Vector3D.Multiply(localDeflection, transformationMatrix);

                    System.Windows.Media.Media3D.Vector3D localRotation = new System.Windows.Media.Media3D.Vector3D(
                        Math.Round(double.Parse(parts[12], CultureInfo.InvariantCulture), 10),
                        Math.Round(double.Parse(parts[13], CultureInfo.InvariantCulture), 10),
                        Math.Round(double.Parse(parts[14], CultureInfo.InvariantCulture), 10)
                    );
                    System.Windows.Media.Media3D.Vector3D globalRotation = System.Windows.Media.Media3D.Vector3D.Multiply(localRotation, transformationMatrix);

                    BEMs bem = new BEMs(time, distanceFromI,
                        Math.Round(double.Parse(parts[3], CultureInfo.InvariantCulture), 10),
                        Math.Round(double.Parse(parts[4], CultureInfo.InvariantCulture), 10),
                        Math.Round(double.Parse(parts[5], CultureInfo.InvariantCulture), 10),
                        Math.Round(double.Parse(parts[6], CultureInfo.InvariantCulture), 10),
                        Math.Round(double.Parse(parts[7], CultureInfo.InvariantCulture), 10),
                        Math.Round(double.Parse(parts[8], CultureInfo.InvariantCulture), 10),
                        globalDeflection.X, globalDeflection.Y, globalDeflection.Z,
                        globalRotation.X, globalRotation.Y, globalRotation.Z
                    );
                    beam.BEMs.Add(bem);
                }
            }
        }
    }
}

Then in my VB.Net desktop app, I attempt to draw the beam extrusions as follow:

Function ExtrudeElements(ByVal bm As oofem_projlib.Beam3D) As Mesh

    Dim startPoint = projoofem.FindNodeById(bm.Node1).ToPoint3D()
    Dim endPoint = projoofem.FindNodeById(bm.Node2).ToPoint3D()

    ' Use the stored local X-axis direction vector
    Dim beamXDirection As Vector3D = bm.LocalXDirection
    beamXDirection.Normalize()

    ' Get the profile points
    Dim profilePoints = bm.CrossSection.ToPoint3D(0)

    ' Create the LinearPath and Region using the profile points
    Dim profileLp As New LinearPath(profilePoints)
    Dim profileRegion As New Region(profileLp)


    ' Determine the Rotation to Align the Profile's Z-axis with the Beam's Direction
    Dim rotationAxis As Vector3D = CrossProduct(Vector3D.AxisZ, beamXDirection)
    Dim angleInRadians As Double = Math.Acos(DotProduct(Vector3D.AxisZ, beamXDirection))

    ' Create the Rotation Transformation
    Dim rotationTrans As New Transformation()
    rotationTrans.Rotation(angleInRadians, rotationAxis)
    profileRegion.TransformBy(rotationTrans)

    ' Translate the profile to the start point of the beam
    Dim translation As New Transformation()
    translation.Translation(startPoint.X, startPoint.Y, startPoint.Z)
    profileRegion.TransformBy(translation)

    ' Use the stored local X-axis for the Extrusion Direction
    Dim extrusionDirection As Vector3D = beamXDirection
    extrusionDirection.Normalize()

    ' Determine the Extrusion Magnitude
    Dim extrusionMagnitude As Double = startPoint.DistanceTo(endPoint)

    ' Extrude the Profile
    Dim extrudedMesh = profileRegion.ExtrudeAsMesh(extrusionMagnitude * extrusionDirection, 0.1, Mesh.natureType.MulticolorSmooth)
    
    Return extrudedMesh

End Function

This works for most part, but some beams are rotated 90 degrees about the local X axis for some reason.

Then, when I attempt to draw the deflected shapes (given as 5 sets of nodes for each beam – start point + 3 internal + end point) things go a bit wrong. The code snippet for this is:

For Each ele In projoofem.Elements

        If ele.ElementType = oofem_projlib.ElementType.Beam3D Then

            Dim bm As oofem_projlib.Beam3D = ele

            Dim bmP1 As Point3D = projoofem.FindNodeById(bm.Node1).ToPoint3D()
            Dim bmP2 As Point3D = projoofem.FindNodeById(bm.Node2).ToPoint3D()

            Dim lstbmPoints As New List(Of Point3D)
            
            For i = 0 To projoofem.nrBeamSegments
                Dim tpt = GetIntermediatePoint(bmP1, bmP2, bm.BEMs(i).DistanceFromI)

                ' Use the global deflection values directly
                Dim xDefl = tpt.X + bm.BEMs(i).D_x * xfrm.tbAmplify.Value
                Dim yDefl = tpt.Y + bm.BEMs(i).D_y * xfrm.tbAmplify.Value
                Dim zDefl = tpt.Z + bm.BEMs(i).D_z * xfrm.tbAmplify.Value

                ''Incorporate rotations(considering small angle approximation)
                'xDefl += bm.BEMs(i).R_z * yDefl - bm.BEMs(i).R_y * zDefl
                'yDefl += bm.BEMs(i).R_x * zDefl - bm.BEMs(i).R_z * xDefl
                'zDefl += bm.BEMs(i).R_y * xDefl - bm.BEMs(i).R_x * yDefl


                lstbmPoints.Add(New Point3D(xDefl, yDefl, zDefl))
                
            Next
            
        End If
    Next

Some deflections are drawn in correctly while others are not. I’m sure the issue is with my coordinate system transformations but I cannot figure out what I’m doing wrong. I’ve spent countless hours trying to solve this. Please help!