diff --git a/VectoCommon/VectoCommon/Utils/DoubleExtensionMethods.cs b/VectoCommon/VectoCommon/Utils/DoubleExtensionMethods.cs index f910671d7f185611ba557de40a23dc012fb82553..10d686da937e2bc3a0cbe76113df7d6e5a78f8f3 100644 --- a/VectoCommon/VectoCommon/Utils/DoubleExtensionMethods.cs +++ b/VectoCommon/VectoCommon/Utils/DoubleExtensionMethods.cs @@ -155,6 +155,18 @@ namespace TUGraz.VectoCommon.Utils return self >= -tolerance; } + /// <summary> + /// Determines whether the specified tolerance is positive within tolerance. + /// </summary> + /// <param name="self">The self.</param> + /// <param name="tolerance">The tolerance.</param> + [DebuggerStepThrough] + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool IsNegative(this double self, double tolerance = Tolerance) + { + return self <= tolerance; + } + /// <summary> /// Checks if a value is between min and max (min <= value <= max) /// </summary> @@ -242,9 +254,9 @@ namespace TUGraz.VectoCommon.Utils return SIBase<Scalar>.Create(value); } - public static SI SI(this double value, UnitInstance si) - { - return new SI(si, value); + public static SI SI(this double value, UnitInstance si) + { + return new SI(si, value); } /// <summary> @@ -279,17 +291,17 @@ namespace TUGraz.VectoCommon.Utils return self.ToString("F" + decimals.Value, CultureInfo.InvariantCulture); } - public static string ToXMLFormat(this ConvertedSI self, uint? decimals = null) - { - decimals = decimals ?? 2; - return ((double)self).ToString("F" + decimals.Value, CultureInfo.InvariantCulture); - } + public static string ToXMLFormat(this ConvertedSI self, uint? decimals = null) + { + decimals = decimals ?? 2; + return ((double)self).ToString("F" + decimals.Value, CultureInfo.InvariantCulture); + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static string ToMinSignificantDigits(this ConvertedSI self, uint? significant = null, uint? decimals = null) - { - return ToMinSignificantDigits((double)self, significant, decimals); - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static string ToMinSignificantDigits(this ConvertedSI self, uint? significant = null, uint? decimals = null) + { + return ToMinSignificantDigits((double)self, significant, decimals); + } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static string ToMinSignificantDigits(this double self, uint? significant = null, uint? decimals = null) @@ -304,6 +316,9 @@ namespace TUGraz.VectoCommon.Utils return self.ToString("F" + Math.Max(significant.Value - scale, decimals.Value), CultureInfo.InvariantCulture); } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool IsNaN(this double self) => double.IsNaN(self); } public static class FloatExtensionMethods diff --git a/VectoCommon/VectoCommon/Utils/VectoMath.cs b/VectoCommon/VectoCommon/Utils/VectoMath.cs index 0932091c61b32ddc024de62493dd84db31879234..e92d493628cb3a156997fd6c3218959a6ed02c47 100644 --- a/VectoCommon/VectoCommon/Utils/VectoMath.cs +++ b/VectoCommon/VectoCommon/Utils/VectoMath.cs @@ -711,22 +711,15 @@ namespace TUGraz.VectoCommon.Utils public bool ContainsInCircumcircle(Point p) { - var p0X = P1.X - p.X; - var p0Y = P1.Y - p.Y; - var p1X = P2.X - p.X; - var p1Y = P2.Y - p.Y; - var p2X = P3.X - p.X; - var p2Y = P3.Y - p.Y; - - var p0Square = p0X * p0X + p0Y * p0Y; - var p1Square = p1X * p1X + p1Y * p1Y; - var p2Square = p2X * p2X + p2Y * p2Y; - - var det01 = p0X * p1Y - p1X * p0Y; - var det12 = p1X * p2Y - p2X * p1Y; - var det20 = p2X * p0Y - p0X * p2Y; - - var result = p0Square * det12 + p1Square * det20 + p2Square * det01; + var ax = P1.X - p.X; + var ay = P1.Y - p.Y; + var bx = P2.X - p.X; + var by = P2.Y - p.Y; + var cx = P3.X - p.X; + var cy = P3.Y - p.Y; + var result = (ax * ax + ay * ay) * (bx * cy - cx * by) + - (bx * bx + by * by) * (ax * cy - cx * ay) + + (cx * cx + cy * cy) * (ax * by - bx * ay); return result > 0; } diff --git a/VectoCore/VectoCore/Models/SimulationComponent/Data/ElectricMotor/EfficiencyMap.cs b/VectoCore/VectoCore/Models/SimulationComponent/Data/ElectricMotor/EfficiencyMap.cs index 6499786445b8a9f7fba85368b221606a3d4227ac..e97a1b0627d3989b9afbe7c62657dffd14e817a4 100644 --- a/VectoCore/VectoCore/Models/SimulationComponent/Data/ElectricMotor/EfficiencyMap.cs +++ b/VectoCore/VectoCore/Models/SimulationComponent/Data/ElectricMotor/EfficiencyMap.cs @@ -25,9 +25,9 @@ namespace TUGraz.VectoCore.Models.SimulationComponent.Data.ElectricMotor var result = new EfficiencyResult(); result.Torque = torque; var value = _efficiencyMapMech2El.Interpolate(torque, angularSpeed); - if (value.HasValue) + if (!value.IsNaN()) { - result.ElectricalPower = value.Value.SI<Watt>(); + result.ElectricalPower = value.SI<Watt>(); return result; } if (allowExtrapolation) diff --git a/VectoCore/VectoCore/Models/SimulationComponent/Data/Engine/FuelConsumptionMap.cs b/VectoCore/VectoCore/Models/SimulationComponent/Data/Engine/FuelConsumptionMap.cs index de041b5ca43812ea56a1a5a21696ff0ee448c798..560d5b7bad656f46c3bd746e60e01bc74dbfafb1 100644 --- a/VectoCore/VectoCore/Models/SimulationComponent/Data/Engine/FuelConsumptionMap.cs +++ b/VectoCore/VectoCore/Models/SimulationComponent/Data/Engine/FuelConsumptionMap.cs @@ -60,8 +60,8 @@ namespace TUGraz.VectoCore.Models.SimulationComponent.Data.Engine var result = new FuelConsumptionResult(); var value = _fuelMap.Interpolate(torque, angularVelocity); - if (value.HasValue) { - result.Value = value.Value.SI(Unit.SI.Kilo.Gramm.Per.Second).Cast<KilogramPerSecond>(); + if (!value.IsNaN()) { + result.Value = value.SI(Unit.SI.Kilo.Gramm.Per.Second).Cast<KilogramPerSecond>(); return result; } diff --git a/VectoCore/VectoCore/Models/SimulationComponent/Data/Engine/WHRPowerMap.cs b/VectoCore/VectoCore/Models/SimulationComponent/Data/Engine/WHRPowerMap.cs index f932d239d8ba0a2bd1525fc6910c03677579f8f5..d3a35bf7a33a2132656e8624ab83285b79228930 100644 --- a/VectoCore/VectoCore/Models/SimulationComponent/Data/Engine/WHRPowerMap.cs +++ b/VectoCore/VectoCore/Models/SimulationComponent/Data/Engine/WHRPowerMap.cs @@ -36,8 +36,8 @@ namespace TUGraz.VectoCore.InputData.Reader.ComponentData var result = new WHRPowerResult(); // delaunay map needs is initialised with rpm, therefore the angularVelocity has to be converted. var value = WHRMap.Interpolate(torque, engineSpeed); - if (value.HasValue) { - result.GeneratedPower = value.Value.SI<Watt>(); + if (!value.IsNaN()) { + result.GeneratedPower = value.SI<Watt>(); return result; } diff --git a/VectoCore/VectoCore/Models/SimulationComponent/Data/Gearbox/TransmissionLossMap.cs b/VectoCore/VectoCore/Models/SimulationComponent/Data/Gearbox/TransmissionLossMap.cs index de9e1b6fb3dcdb50ba3e6532964fa132e80dbae3..9a5d977006d8ac60ad9d4eb69ff89ac2ea3eaa34 100644 --- a/VectoCore/VectoCore/Models/SimulationComponent/Data/Gearbox/TransmissionLossMap.cs +++ b/VectoCore/VectoCore/Models/SimulationComponent/Data/Gearbox/TransmissionLossMap.cs @@ -93,12 +93,12 @@ namespace TUGraz.VectoCore.Models.SimulationComponent.Data.Gearbox var result = new LossMapResult(); var torqueLoss = _lossMap.Interpolate(outAngularVelocity.Value() * _ratio, outTorque.Value() / _ratio); - if (!torqueLoss.HasValue) { + if (torqueLoss.IsNaN()) { torqueLoss = _lossMap.Extrapolate(outAngularVelocity.Value() * _ratio, outTorque.Value() / _ratio); result.Extrapolated = true; } - result.Value = torqueLoss.Value.SI<NewtonMeter>(); + result.Value = torqueLoss.SI<NewtonMeter>(); Log.Debug("GearboxLoss {0}: {1}, outAngularVelocity: {2}, outTorque: {3}", GearName, torqueLoss, outAngularVelocity, outTorque); @@ -122,13 +122,13 @@ namespace TUGraz.VectoCore.Models.SimulationComponent.Data.Gearbox public NewtonMeter GetOutTorque(PerSecond inAngularVelocity, NewtonMeter inTorque, bool allowExtrapolation = false) { var torqueLoss = _invertedLossMap.Interpolate(inAngularVelocity.Value(), inTorque.Value()); - if (torqueLoss.HasValue) { - return (inTorque - torqueLoss.Value.SI<NewtonMeter>()) * _ratio; + if (!torqueLoss.IsNaN()) { + return (inTorque - torqueLoss.SI<NewtonMeter>()) * _ratio; } if (allowExtrapolation) { torqueLoss = _invertedLossMap.Extrapolate(inAngularVelocity.Value(), inTorque.Value()); - return (inTorque - torqueLoss.Value.SI<NewtonMeter>()) * _ratio; + return (inTorque - torqueLoss.SI<NewtonMeter>()) * _ratio; } throw new VectoException("TransmissionLossMap {0}: Interpolation failed. inTorque: {1}, inAngularVelocity: {2}", diff --git a/VectoCore/VectoCore/Utils/DelaunayMap.cs b/VectoCore/VectoCore/Utils/DelaunayMap.cs index c9e483cd3203f182657487dcd11f57b667618fb8..86b43e02e0ac81e7eb0967b3b4ce8596251738a5 100644 --- a/VectoCore/VectoCore/Utils/DelaunayMap.cs +++ b/VectoCore/VectoCore/Utils/DelaunayMap.cs @@ -37,7 +37,6 @@ using System.IO; using System.Linq; using System.Runtime.CompilerServices; using System.Windows.Forms.DataVisualization.Charting; -using Newtonsoft.Json; using TUGraz.VectoCommon.Exceptions; using TUGraz.VectoCommon.Models; using TUGraz.VectoCommon.Utils; @@ -71,8 +70,7 @@ namespace TUGraz.VectoCore.Utils public IReadOnlyCollection<Point> Entries { - get - { + get { var retVal = new Point[_points.Count]; var i = 0; foreach (var pt in _points) { @@ -103,43 +101,42 @@ namespace TUGraz.VectoCore.Utils _maxY = _points.Max(p => p.Y); _minX = _points.Min(p => p.X); _minY = _points.Min(p => p.Y); - _points = - _points.Select(p => new Point((p.X - _minX) / (_maxX - _minX), (p.Y - _minY) / (_maxY - _minY), p.Z)).ToList(); + _points = _points.Select(p => new Point((p.X - _minX) / (_maxX - _minX), (p.Y - _minY) / (_maxY - _minY), p.Z)).ToArray(); var superTriangle = new Triangle(new Point(-1, -1), new Point(4, -1), new Point(-1, 4)); var triangles = new List<Triangle> { superTriangle }; - var pointCount = 0; - - var points = _points.ToArray(); - - // iteratively add each point into the correct triangle and split up the triangle - foreach (var point in points) { + var n = 3; //starts with 3 points on supertriangle + // iteratively add each point into the correct triangle and split up the triangle + var uniqueEdges = new HashSet<Edge>(); + foreach (var point in _points) { // If the vertex lies inside the circumcircle of a triangle, the edges of this triangle are // added to the edge buffer and the triangle is removed from list. // Remove duplicate edges. This leaves the convex hull of the edges. // The edges in this convex hull are oriented counterclockwise! - var newTriangles = triangles.Select((t, i) => Tuple.Create(i, t, t.ContainsInCircumcircle(point))) - .Where(t => t.Item3) - .Reverse() - .SelectMany(t => { - triangles.RemoveAt(t.Item1); - return t.Item2.GetEdges(); - }) - .GroupBy(edge => edge) - .Where(group => group.Count() == 1) - .Select(group => new Triangle(group.Key.P1, group.Key.P2, point)).ToList(); - - triangles.AddRange(newTriangles); - - //DrawGraph(pointCount, triangles, superTriangle, xmin, xmax, ymin, ymax, point); - pointCount++; + uniqueEdges.Clear(); + for (var i = triangles.Count - 1; i >= 0; i--) { + if (triangles[i].ContainsInCircumcircle(point)) { + foreach (var edge in triangles[i].GetEdges()) { + if (uniqueEdges.Contains(edge)) { + uniqueEdges.Remove(edge); + } else { + uniqueEdges.Add(edge); + } + } + triangles.RemoveAt(i); + } + } + foreach (var edge in uniqueEdges) { + triangles.Add(new Triangle(point, edge.P1, edge.P2)); + } + n++; - // check invariant: m = 2n-2-k + // triangle invariant: m=2n-2-k // m...triangle count - // n...point count (pointCount +3 points on the supertriangle) + // n...point count (pointCount including 3 points on the supertriangle) // k...points on convex hull (exactly 3 --> supertriangle) - if (triangles.Count != 2 * (pointCount + 3) - 2 - 3) { + if (triangles.Count != 2 * n - 2 - 3) { throw new VectoException( "{0} Delaunay-Triangulation invariant violated! Triangle count and point count doesn't fit together.", _mapName); } @@ -229,10 +226,8 @@ namespace TUGraz.VectoCore.Utils } } - public double? Interpolate(SI x, SI y) - { - return Interpolate(x.Value(), y.Value()); - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public double Interpolate(SI x, SI y) => Interpolate(x.Value(), y.Value()); /// <summary> /// Interpolates the value of an point in the delaunay map. @@ -242,7 +237,7 @@ namespace TUGraz.VectoCore.Utils /// <returns>a value if interpolation is successfull, /// null if interpolation has failed.</returns> [MethodImpl(MethodImplOptions.AggressiveInlining)] - public double? Interpolate(double x, double y) + public double Interpolate(double x, double y) { if (_triangles == null) { throw new VectoException("Interpolation not possible. Call DelaunayMap.Triangulate first."); @@ -250,7 +245,7 @@ namespace TUGraz.VectoCore.Utils x = (x - _minX) / (_maxX - _minX); y = (y - _minY) / (_maxY - _minY); - + var i = 0; while (i < _triangles.Length && !_triangles[i].IsInside(x, y, true)) { i++; @@ -260,15 +255,27 @@ namespace TUGraz.VectoCore.Utils while (i < _triangles.Length && !_triangles[i].IsInside(x, y, false)) { i++; } - } - if (i == _triangles.Length) { - return null; + if (i == _triangles.Length) { + return double.NaN; + } } var tr = _triangles[i]; - var plane = new Plane(tr); - return (plane.W - plane.X * x - plane.Y * y) / plane.Z; + var abX = tr.P2.X - tr.P1.X; + var abY = tr.P2.Y - tr.P1.Y; + var abZ = tr.P2.Z - tr.P1.Z; + + var acX = tr.P3.X - tr.P1.X; + var acY = tr.P3.Y - tr.P1.Y; + var acZ = tr.P3.Z - tr.P1.Z; + + var X = abY * acZ - abZ * acY; + var Y = abZ * acX - abX * acZ; + var Z = abX * acY - abY * acX; + var W = tr.P1.X * X + tr.P1.Y * Y + tr.P1.Z * Z; + + return (W - X * x - Y * y) / Z; } public double Extrapolate(SI x, SI y) diff --git a/VectoCore/VectoCoreTest/Algorithms/DelaunayMapTest.cs b/VectoCore/VectoCoreTest/Algorithms/DelaunayMapTest.cs index e6017218e5fa8c835f64c2b1fed8c83ee0bacf1b..baae14e9217eca3e6a489ea70b347eb47a03afbd 100644 --- a/VectoCore/VectoCoreTest/Algorithms/DelaunayMapTest.cs +++ b/VectoCore/VectoCoreTest/Algorithms/DelaunayMapTest.cs @@ -94,10 +94,10 @@ namespace TUGraz.VectoCore.Tests.Algorithms AssertHelper.AreRelativeEqual(1.5, map.Interpolate(0, 0.75)); // extrapolation (should fail) - Assert.IsNull(map.Interpolate(1, 1)); - Assert.IsNull(map.Interpolate(-1, -1)); - Assert.IsNull(map.Interpolate(1, -1)); - Assert.IsNull(map.Interpolate(-1, 1)); + Assert.IsNaN(map.Interpolate(1, 1)); + Assert.IsNaN(map.Interpolate(-1, -1)); + Assert.IsNaN(map.Interpolate(1, -1)); + Assert.IsNaN(map.Interpolate(-1, 1)); } public void Test_DelaunayMapPlane()