Fixed computation of aspect ration

I backported the code from cpacs creator,
which already solved the problem.

Fixes #827

src/api/tigl.h

@@ -146,6 +146,20 @@ enum TiglBoolean
 typedef enum TiglBoolean TiglBoolean;
 
 
+/**
+ \ingroup Enums
+  Definition of Axis used in TIGL.
+*/
+enum TiglAxis
+{
+    TIGL_NO_AXIS = 0,
+    TIGL_X_AXIS = 1,
+    TIGL_Y_AXIS = 2,
+    TIGL_Z_AXIS = 3
+};
+
+typedef enum TiglAxis TiglAxis;
+
 /**
  \ingroup Enums
   Definition of Symmetry Axis used in TIGL.

src/math/tiglmathfunctions.cpp

@@ -452,4 +452,9 @@ double Interpolate(const std::vector<double>& xdata, const std::vector<double>&
     return y;
 }
 
+bool isNear(double a, double b, double epsilon)
+{
+    return (fabs(a - b) <= epsilon);
+}
+
 } // namespace tigl

src/math/tiglmathfunctions.h

@@ -23,6 +23,7 @@
 #include "tigl_internal.h"
 #include <vector>
 #include <math_Vector.hxx>
+#include <Precision.hxx>
 #include "tiglMatrix.h"
 
 namespace tigl 
@@ -127,6 +128,14 @@ TIGL_EXPORT double cstcurve(const double& N1, const double& N2, const std::vecto
  */
 TIGL_EXPORT double cstcurve_deriv(const double& N1, const double& N2, const std::vector<double>& B, const double& T, const int& n, const double& x);
 
+/**
+ * Return true if the value of a is similar to b
+ * @param a
+ * @param b
+ * @param epsilon
+ * @return
+ */
+TIGL_EXPORT bool isNear(double a, double b, double epsilon = Precision::Confusion());
 
 /**
  * 1D Function interface accepting one parameter t and returning

src/wing/CCPACSWing.cpp

@@ -41,6 +41,7 @@
 #include "CTiglBSplineAlgorithms.h"
 #include "CCPACSControlSurfaces.h"
 #include "CCPACSTrailingEdgeDevice.h"
+#include "TiglWingHelperFunctions.h"
 
 #include "BRepOffsetAPI_ThruSections.hxx"
 #include "BRepAlgoAPI_Fuse.hxx"
@@ -68,16 +69,6 @@ namespace tigl
 
 namespace
 {
-    inline double max(double a, double b)
-    {
-        return a > b? a : b;
-    }
-
-    inline double min(double a, double b)
-    {
-        return a < b? a : b;
-    }
-
     // Returns the index of the maximum value
     int maxIndex(double x, double y, double z)
     {
@@ -547,6 +538,35 @@ double CCPACSWing::GetReferenceArea(TiglSymmetryAxis symPlane) const
     return refArea;
 }
 
+double CCPACSWing::GetReferenceArea() const
+{
+    TiglAxis spanDir = winghelper::GetMajorDirection(*this);
+    TiglAxis deepDir = winghelper::GetDeepDirection(*this);
+
+    if (spanDir == TIGL_Y_AXIS && deepDir == TIGL_X_AXIS) {
+        return GetReferenceArea(TIGL_X_Y_PLANE);
+    }
+    else if (spanDir == TIGL_Y_AXIS && deepDir == TIGL_Z_AXIS) {
+        return GetReferenceArea(TIGL_Y_Z_PLANE);
+    }
+    else if (spanDir == TIGL_Z_AXIS && deepDir == TIGL_X_AXIS) {
+        return GetReferenceArea(TIGL_X_Z_PLANE);
+    }
+    else if (spanDir == TIGL_Z_AXIS && deepDir == TIGL_Y_AXIS) {
+        return GetReferenceArea(TIGL_Y_Z_PLANE);
+    }
+    else if (spanDir == TIGL_X_AXIS && deepDir == TIGL_Z_AXIS) {
+        return GetReferenceArea(TIGL_X_Z_PLANE);
+    }
+    else if (spanDir == TIGL_X_AXIS && deepDir == TIGL_Y_AXIS) {
+        return GetReferenceArea(TIGL_X_Y_PLANE);
+    }
+    else {
+       LOG(ERROR) << "CCPACSWing::GetReferenceArea: Unexpected pair of major direction and deep direction.";
+       return 0.0;
+    }
+}
+
 
 double CCPACSWing::GetWettedArea(TopoDS_Shape parent) const
 {
@@ -663,7 +683,21 @@ double CCPACSWing::GetWingspan() const
 //     s: half span; A_half: Reference area of wing without symmetrical wing
 double CCPACSWing::GetAspectRatio() const
 {
-    return 2.0*(pow_int(GetWingspan(),2)/GetReferenceArea(GetSymmetryAxis()));
+    auto refArea = GetReferenceArea();
+
+    if ( isNear(refArea, 0.) ) {
+        LOG(WARNING) << "Wing area is close to zero, thus the AR is not computed and 0 is returned.";
+        return 0;
+    }
+
+    auto halfSpan = GetWingspan();
+
+    // If the wing has symmetry defined, the full span was returned
+    if (GetSymmetryAxis() != TIGL_NO_SYMMETRY) {
+        halfSpan *= 0.5;
+    }
+
+    return 2.0*halfSpan*halfSpan/refArea;
 }
 
 /**

src/wing/CCPACSWing.h

@@ -125,6 +125,9 @@ friend class CTiglWingBuilder;
     // of each wing segment by projecting the wing segments into the plane defined by the user
     TIGL_EXPORT double GetReferenceArea(TiglSymmetryAxis symPlane) const;
 
+    // Returns the reference area of the wing in the plane normal to the major direction
+    TIGL_EXPORT double GetReferenceArea() const;
+
     // Returns wetted Area
     TIGL_EXPORT double GetWettedArea(TopoDS_Shape parent) const;
 

src/wing/TiglWingHelperFunctions.cpp

@@ -0,0 +1,118 @@
+#include "TiglWingHelperFunctions.h"
+
+#include "tigl.h"
+#include "CCPACSWingSegment.h"
+
+namespace tigl
+{
+
+namespace winghelper
+{
+
+bool HasShape(const tigl::CCPACSWing& wing)
+{
+    return wing.GetSegmentCount() > 0;
+
+}
+
+TiglAxis GetDeepDirection(const tigl::CCPACSWing& wing)
+{
+    if (!HasShape(wing)) {
+        LOG(WARNING) << "CTiglWingHelper::GetDeepDirection: This wing has no shape -> impossible to determine the "
+                        "direction properly. The default direction will be returned";
+        return TIGL_X_AXIS;
+    }
+
+    gp_XYZ cumulatedDepthDirection(0, 0, 0);
+    for (int i = 1; i <= wing.GetSegmentCount(); ++i) {
+        const CCPACSWingSegment& segment = wing.GetSegment(i);
+        gp_XYZ dirDepth                  = segment.GetChordPoint(0, 1).XYZ() - segment.GetChordPoint(0, 0).XYZ();
+        dirDepth                         = gp_XYZ(fabs(dirDepth.X()), fabs(dirDepth.Y()), fabs(dirDepth.Z()));
+        cumulatedDepthDirection += dirDepth;
+    }
+    const CCPACSWingSegment& outerSegment = wing.GetSegment(wing.GetSegmentCount());
+    gp_XYZ dirDepth = outerSegment.GetChordPoint(1, 1).XYZ() - outerSegment.GetChordPoint(1, 0).XYZ();
+    dirDepth        = gp_XYZ(fabs(dirDepth.X()), fabs(dirDepth.Y()), fabs(dirDepth.Z()));
+    cumulatedDepthDirection += dirDepth;
+
+    switch (GetMajorDirection(wing)) {
+    case TIGL_Y_AXIS:
+        return cumulatedDepthDirection.X() >= cumulatedDepthDirection.Z() ? TiglAxis::TIGL_X_AXIS
+                                                                          : TiglAxis::TIGL_Z_AXIS;
+    case TIGL_Z_AXIS:
+        return cumulatedDepthDirection.X() >= cumulatedDepthDirection.Y() ? TiglAxis::TIGL_X_AXIS
+                                                                          : TiglAxis::TIGL_Y_AXIS;
+    case TIGL_X_AXIS:
+        return cumulatedDepthDirection.Z() >= cumulatedDepthDirection.Y() ? TiglAxis::TIGL_Z_AXIS
+                                                                          : TiglAxis::TIGL_Y_AXIS;
+    default:
+        return TiglAxis ::TIGL_X_AXIS;
+    }
+}
+
+TiglAxis GetMajorDirection(const CCPACSWing& wing)
+{
+    if (!HasShape(wing)) {
+        LOG(WARNING) << "CTiglWingHelper::GetMajorDirection: This wing has no shape -> impossible to determine the "
+                        "direction properly. The default direction will be returned";
+        return TIGL_Y_AXIS;
+    }
+
+    switch (wing.GetSymmetryAxis()) {
+    case TIGL_X_Y_PLANE:
+        return TiglAxis::TIGL_Z_AXIS;
+    case TIGL_X_Z_PLANE:
+        return TiglAxis::TIGL_Y_AXIS;
+    case TIGL_Y_Z_PLANE:
+        return TiglAxis ::TIGL_X_AXIS;
+    default:
+        // heuristic to find the best major axis
+        // first find the deep axis , then chose the major axis between the two left axis
+        gp_XYZ cumulatedSpanDirection(0, 0, 0);
+        gp_XYZ cumulatedDepthDirection(0, 0, 0);
+        for (int i = 1; i <= wing.GetSegmentCount(); ++i) {
+            const CCPACSWingSegment& segment = wing.GetSegment(i);
+            gp_XYZ dirSpan                   = segment.GetChordPoint(1, 0).XYZ() - segment.GetChordPoint(0, 0).XYZ();
+            gp_XYZ dirDepth                  = segment.GetChordPoint(0, 1).XYZ() - segment.GetChordPoint(0, 0).XYZ();
+            dirSpan  = gp_XYZ(fabs(dirSpan.X()), fabs(dirSpan.Y()), fabs(dirSpan.Z())); // why we use abs value?
+            dirDepth = gp_XYZ(fabs(dirDepth.X()), fabs(dirDepth.Y()), fabs(dirDepth.Z()));
+            cumulatedSpanDirection += dirSpan;
+            cumulatedDepthDirection += dirDepth;
+        }
+        const CCPACSWingSegment& outerSegment = wing.GetSegment(wing.GetSegmentCount());
+        gp_XYZ dirDepth = outerSegment.GetChordPoint(1, 1).XYZ() - outerSegment.GetChordPoint(1, 0).XYZ();
+        dirDepth        = gp_XYZ(fabs(dirDepth.X()), fabs(dirDepth.Y()), fabs(dirDepth.Z()));
+        cumulatedDepthDirection += dirDepth;
+
+        int depthIndex = 0;
+        if (cumulatedDepthDirection.X() >= cumulatedDepthDirection.Y() &&
+            cumulatedDepthDirection.X() >= cumulatedDepthDirection.Z()) {
+            depthIndex = 0;
+        }
+        else if (cumulatedDepthDirection.Y() >= cumulatedDepthDirection.X() &&
+                 cumulatedDepthDirection.Y() >= cumulatedDepthDirection.Z()) {
+            depthIndex = 1;
+        }
+        else {
+            depthIndex = 2;
+        }
+
+        switch (depthIndex) {
+        case 0:
+            return cumulatedSpanDirection.Y() >= cumulatedSpanDirection.Z() ? TiglAxis::TIGL_Y_AXIS
+                                                                            : TiglAxis::TIGL_Z_AXIS;
+        case 1:
+            return cumulatedSpanDirection.X() >= cumulatedSpanDirection.Z() ? TiglAxis::TIGL_X_AXIS
+                                                                            : TiglAxis::TIGL_Z_AXIS;
+        case 2:
+            return cumulatedSpanDirection.X() >= cumulatedSpanDirection.Y() ? TiglAxis::TIGL_X_AXIS
+                                                                            : TiglAxis::TIGL_Y_AXIS;
+        default:
+            return TiglAxis ::TIGL_Y_AXIS;
+        }
+    }
+}
+
+} // namespace winghelper
+
+} // namespace tigl

src/wing/TiglWingHelperFunctions.h

@@ -0,0 +1,24 @@
+#ifndef TIGLWINGHELPERFUNCTIONS_H
+#define TIGLWINGHELPERFUNCTIONS_H
+
+#include <CCPACSWing.h>
+
+namespace tigl
+{
+
+namespace winghelper
+{
+
+ // Returns the deep direction of the wing
+TiglAxis GetDeepDirection(const tigl::CCPACSWing& wing);
+
+// Returns the major direction of the wing (correspond to the span direction)
+// @Details: If a symmetry plan is set, the major direction is normal to the symmetry plan,
+// otherwise, an heuristic is used to find out the best span axis candidate.
+TiglAxis GetMajorDirection(const tigl::CCPACSWing& wing);
+
+} // namespace winghelper
+
+} // namespace tigl
+
+#endif // TIGLWINGHELPERFUNCTIONS_H

tests/unittests/tiglWing.cpp

@@ -298,16 +298,17 @@ TEST_F(TiglWing, tiglWingGetSpanVTP)
     ASSERT_GT(span, 5.8);
 }
 
+/// Regression reported by issue #827
 TEST_F(TiglWing, tiglGetAspectRatio)
 {
     tigl::CCPACSConfigurationManager& manager = tigl::CCPACSConfigurationManager::GetInstance();
     tigl::CCPACSConfiguration& config = manager.GetConfiguration(tiglHandle);
 
     const auto& wing = config.GetWing(1);
-    auto ar = wing.GetAspectRatio();
-    EXPECT_NEAR(9.4, ar, 0.1);
-}
 
+    // The expected value was provided by the issue report
+    EXPECT_NEAR(9.4031, wing.GetAspectRatio(), 1e-4);
+}
 
 TEST_F(WingSimple, wingGetMAC_success)
 {