From 559abbbfc2a5c233c8264aa515cea1136fb8c0cd Mon Sep 17 00:00:00 2001 From: Eduard Date: Sun, 13 Sep 2020 17:42:06 +0200 Subject: [PATCH] remove math.net from main project --- NLES.Tests/Builder.cs | 70 ++++++++-------- NLES.Tests/LinearSolverForTesting.cs | 26 ++++-- NLES.Tests/NLES.Tests.csproj | 2 +- NLES.Tests/TestSolverResult.cs | 58 +++++++------ NLES/Builder/SolverBuilder.cs | 8 +- NLES/Builder/SolverBuilderIncrementalLoad.cs | 7 +- NLES/Builder/SolverBuilderStructure.cs | 18 ++-- NLES/Contracts/LoadState.cs | 7 +- NLES/Correction/CorrectionSchemeArcLength.cs | 9 +- NLES/Correction/CorrectionSchemeStandard.cs | 8 +- NLES/Correction/Corrector.cs | 16 ++-- NLES/Correction/ICorrectionScheme.cs | 7 +- NLES/Correction/Methods/RestoringMethod.cs | 9 +- NLES/ILinearSolver.cs | 5 +- NLES/LoadIncrementalState.cs | 7 +- NLES/NLES.csproj | 1 - NLES/NonLinearSolver.cs | 4 +- NLES/Prediction/IPredictionScheme.cs | 5 +- NLES/Prediction/PredictionSchemeArcLength.cs | 3 +- NLES/Prediction/PredictionSchemeStandard.cs | 5 +- NLES/Prediction/Predictor.cs | 13 ++- NLES/StructureInfo.cs | 10 +-- NLES/Vector.cs | 88 ++++++++++++++++++++ 23 files changed, 230 insertions(+), 156 deletions(-) create mode 100644 NLES/Vector.cs diff --git a/NLES.Tests/Builder.cs b/NLES.Tests/Builder.cs index 2ae9db8..beb0557 100644 --- a/NLES.Tests/Builder.cs +++ b/NLES.Tests/Builder.cs @@ -1,15 +1,13 @@ using System; using System.Linq; - -using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearAlgebra.Double; - using NLES; using NLES.Correction; using NLES.Correction.Methods; using NLES.Tests; using Xunit; +using Vector = NLES.Vector; namespace NonLinearEquationsSolver { @@ -19,14 +17,14 @@ public class Builder public void RestoringArcLengthSolverIsCorrectlyBuild() { //Given - static Vector Reaction(Vector u) => new DenseVector(2); + static Vector Reaction(Vector u) => new Vector(2); - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2)); - Vector initialLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; - Vector initialDisplacement = new DenseVector(2); + Vector initialLoad = new Vector(2) { [0] = 1, [1] = 1 }; + Vector initialDisplacement = new Vector(2); double initialLoadFactor = 1; - Vector referenceLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; + Vector referenceLoad = new Vector(2) { [0] = 1, [1] = 1 }; double dispTol = 1e-3; double eqTol = 1e-3; double enTol = 1e-3; @@ -62,14 +60,14 @@ static ILinearSolver Stiffness(Vector u) public void AngleArcLengthSolverIsCorrectlyBuild() { //Given - static Vector Reaction(Vector u) => new DenseVector(2); + static Vector Reaction(Vector u) => new Vector(2); - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2)); - Vector initialLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; - Vector initialDisplacement = new DenseVector(2); + Vector initialLoad = new Vector(2) { [0] = 1, [1] = 1 }; + Vector initialDisplacement = new Vector(2); double initialLoadFactor = 1; - Vector referenceLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; + Vector referenceLoad = new Vector(2) { [0] = 1, [1] = 1 }; double dispTol = 1e-3; double eqTol = 1e-3; double enTol = 1e-3; @@ -105,14 +103,14 @@ static ILinearSolver Stiffness(Vector u) public void StandardSolverIsCorrectlyBuild() { //Given - static Vector Reaction(Vector u) => new DenseVector(2); + static Vector Reaction(Vector u) => new Vector(2); - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2)); - Vector initialLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; - Vector initialDisplacement = new DenseVector(2); + Vector initialLoad = new Vector(2) { [0] = 1, [1] = 1 }; + Vector initialDisplacement = new Vector(2); double initialLoadFactor = 1; - Vector referenceLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; + Vector referenceLoad = new Vector(2) { [0] = 1, [1] = 1 }; double dispTol = 1e-3; double eqTol = 1e-3; double enTol = 1e-3; @@ -144,10 +142,10 @@ static ILinearSolver Stiffness(Vector u) public void StandardSolverWithoutReactionDoesNotLaunchException() { //Given - Vector initialLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; - Vector initialDisplacement = new DenseVector(2); + Vector initialLoad = new Vector(2) { [0] = 1, [1] = 1 }; + Vector initialDisplacement = new Vector(2); double initialLoadFactor = 1; - Vector referenceLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; + Vector referenceLoad = new Vector(2) { [0] = 1, [1] = 1 }; double dispTol = 1e-3; double eqTol = 1e-3; double enTol = 1e-3; @@ -177,14 +175,14 @@ public void StandardSolverWithoutReactionDoesNotLaunchException() public void StandardSolverWithoutDefinedStopConditionsHasDefaultStopConditions() { //Given - static Vector Reaction(Vector u) => new DenseVector(2); + static Vector Reaction(Vector u) => new Vector(2); - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2)); - Vector initialLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; - Vector initialDisplacement = new DenseVector(2); + Vector initialLoad = new Vector(2) { [0] = 1, [1] = 1 }; + Vector initialDisplacement = new Vector(2); double initialLoadFactor = 1; - Vector referenceLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; + Vector referenceLoad = new Vector(2) { [0] = 1, [1] = 1 }; double dlambda = 0.1; // When @@ -210,14 +208,14 @@ static ILinearSolver Stiffness(Vector u) public void StandardSolverWithoutDefinedLoadFactorIncrementLaunchesException() { //Given - static Vector Reaction(Vector u) => new DenseVector(2); + static Vector Reaction(Vector u) => new Vector(2); - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2)); - Vector initialLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; - Vector initialDisplacement = new DenseVector(2); + Vector initialLoad = new Vector(2) { [0] = 1, [1] = 1 }; + Vector initialDisplacement = new Vector(2); double initialLoadFactor = 1; - Vector referenceLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; + Vector referenceLoad = new Vector(2) { [0] = 1, [1] = 1 }; bool exceptionLaunched = false; try @@ -241,14 +239,14 @@ static ILinearSolver Stiffness(Vector u) public void ArcLengthSolverWithoutDefinedRadiusLaunchesException() { //Given - static Vector Reaction(Vector u) => new DenseVector(2); + static Vector Reaction(Vector u) => new Vector(2); - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2)); - Vector initialLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; - Vector initialDisplacement = new DenseVector(2); + Vector initialLoad = new Vector(2) { [0] = 1, [1] = 1 }; + Vector initialDisplacement = new Vector(2); double initialLoadFactor = 1; - Vector referenceLoad = new DenseVector(2) { [0] = 1, [1] = 1 }; + Vector referenceLoad = new Vector(2) { [0] = 1, [1] = 1 }; bool exceptionLaunched = false; try diff --git a/NLES.Tests/LinearSolverForTesting.cs b/NLES.Tests/LinearSolverForTesting.cs index d4231a0..cfa0517 100644 --- a/NLES.Tests/LinearSolverForTesting.cs +++ b/NLES.Tests/LinearSolverForTesting.cs @@ -1,7 +1,5 @@ -using System; -using System.Collections.Generic; -using System.Text; -using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; namespace NLES.Tests { @@ -11,7 +9,23 @@ internal class LinearSolverForTesting : ILinearSolver public LinearSolverForTesting(Matrix matrix) => this.matrix = matrix; - public Vector Solve(Vector input) - => matrix.Solve(input); + public Vector Solve(Vector input) + { + Vector v = new DenseVector(input.Dimension); + for (int i = 0; i < input.Dimension; i++) + { + v[i] = input[i]; + } + + Vector sol = matrix.Solve(v); + + Vector solution = new Vector(input.Dimension); + for (int i = 0; i < input.Dimension; i++) + { + solution[i] = sol[i]; + } + + return solution; + } } } diff --git a/NLES.Tests/NLES.Tests.csproj b/NLES.Tests/NLES.Tests.csproj index a6b3190..09d3fe1 100644 --- a/NLES.Tests/NLES.Tests.csproj +++ b/NLES.Tests/NLES.Tests.csproj @@ -5,7 +5,7 @@ - + diff --git a/NLES.Tests/TestSolverResult.cs b/NLES.Tests/TestSolverResult.cs index 8d1fe81..237209f 100644 --- a/NLES.Tests/TestSolverResult.cs +++ b/NLES.Tests/TestSolverResult.cs @@ -1,14 +1,12 @@ using System; using System.Collections.Generic; using System.Linq; - -using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearAlgebra.Double; - using NLES; using NLES.Contracts; using NLES.Tests; using Xunit; +using Vector = NLES.Vector; namespace NonLinearEquationsSolver { @@ -21,11 +19,11 @@ public class TestSolverResults public void SolveLinearFunction() { // ARRANGE - static Vector Reaction(Vector u) => new DenseVector(2) { [0] = u[0], [1] = u[0] + 2 * u[1] }; + static Vector Reaction(Vector u) => new Vector(2) { [0] = u[0], [1] = u[0] + 2 * u[1] }; - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2) { [0, 0] = 1, [1, 0] = 1, [1, 1] = 2 }); - DenseVector force = new DenseVector(2) { [0] = 1, [1] = 3 }; + Vector force = new Vector(2) { [0] = 1, [1] = 3 }; NonLinearSolver Solver = NonLinearSolver.Builder .Solve(2, Reaction, Stiffness) .Under(force) @@ -42,11 +40,11 @@ static ILinearSolver Stiffness(Vector u) public void SolveLinearFunctionSmallIncrements() { // ARRANGE - static Vector Reaction(Vector u) => new DenseVector(2) { [0] = u[0], [1] = u[0] + 2 * u[1] }; + static Vector Reaction(Vector u) => new Vector(2) { [0] = u[0], [1] = u[0] + 2 * u[1] }; - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2) { [0, 0] = 1, [1, 0] = 1, [1, 1] = 2 }); - DenseVector force = new DenseVector(2) { [0] = 1, [1] = 3 }; + Vector force = new Vector(2) { [0] = 1, [1] = 3 }; double inc = 1e-2; NonLinearSolver Solver = NonLinearSolver.Builder .Solve(2, Reaction, Stiffness) @@ -66,11 +64,11 @@ static ILinearSolver Stiffness(Vector u) public void SolveLinearFunctionArcLength() { // ARRANGE - static Vector Reaction(Vector u) => new DenseVector(2) { [0] = u[0], [1] = u[0] + 2 * u[1] }; + static Vector Reaction(Vector u) => new Vector(2) { [0] = u[0], [1] = u[0] + 2 * u[1] }; - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2) { [0, 0] = 1, [1, 0] = 1, [1, 1] = 2 }); - DenseVector force = new DenseVector(2) { [0] = 1, [1] = 3 }; + Vector force = new Vector(2) { [0] = 1, [1] = 3 }; double radius = 1e-2; NonLinearSolver Solver = NonLinearSolver.Builder .Solve(2, Reaction, Stiffness) @@ -89,13 +87,13 @@ static ILinearSolver Stiffness(Vector u) public void SolveQuadraticFunction() { // ARRANGE - static Vector Reaction(Vector u) => new DenseVector(2) + static Vector Reaction(Vector u) => new Vector(2) { [0] = u[0] * u[0] + 2 * u[1] * u[1], [1] = 2 * u[0] * u[0] + u[1] * u[1] }; - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2) { [0, 0] = 2 * u[0], @@ -103,11 +101,11 @@ static ILinearSolver Stiffness(Vector u) [1, 0] = 4 * u[0], [1, 1] = 2 * u[1] }); - DenseVector force = new DenseVector(2) { [0] = 3, [1] = 3 }; + Vector force = new Vector(2) { [0] = 3, [1] = 3 }; NonLinearSolver Solver = NonLinearSolver.Builder .Solve(2, Reaction, Stiffness) .Under(force) - .WithInitialConditions(0.1, DenseVector.Create(2, 0), DenseVector.Create(2, 1)) + .WithInitialConditions(0.1, new Vector(2, 0), new Vector(2, 1)) .Build(); // ACT @@ -121,13 +119,13 @@ static ILinearSolver Stiffness(Vector u) public void SolveQuadraticFunctionSmallIncrements() { // ARRANGE - static Vector Reaction(Vector u) => new DenseVector(2) + static Vector Reaction(Vector u) => new Vector(2) { [0] = u[0] * u[0] + 2 * u[1] * u[1], [1] = 2 * u[0] * u[0] + u[1] * u[1] }; - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2) { [0, 0] = 2 * u[0], @@ -135,11 +133,11 @@ static ILinearSolver Stiffness(Vector u) [1, 0] = 4 * u[0], [1, 1] = 2 * u[1] }); - DenseVector force = new DenseVector(2) { [0] = 3, [1] = 3 }; + Vector force = new Vector(2) { [0] = 3, [1] = 3 }; NonLinearSolver Solver = NonLinearSolver.Builder .Solve(2, Reaction, Stiffness) .Under(force) - .WithInitialConditions(0.1, DenseVector.Create(2, 0), DenseVector.Create(2, 1)) + .WithInitialConditions(0.1, new Vector(2, 0), new Vector(2, 1)) .UsingStandardNewtonRaphsonScheme(0.01) .Build(); @@ -154,13 +152,13 @@ static ILinearSolver Stiffness(Vector u) public void SolveQuadraticFunctionArcLength() { // ARRANGE - static Vector Reaction(Vector u) => new DenseVector(2) + static Vector Reaction(Vector u) => new Vector(2) { [0] = u[0] * u[0] + 2 * u[1] * u[1], [1] = 2 * u[0] * u[0] + u[1] * u[1] }; - static ILinearSolver Stiffness(Vector u) + static ILinearSolver Stiffness(Vector u) => new LinearSolverForTesting(new DenseMatrix(2, 2) { [0, 0] = 2 * u[0], @@ -168,11 +166,11 @@ static ILinearSolver Stiffness(Vector u) [1, 0] = 4 * u[0], [1, 1] = 2 * u[1] }); - DenseVector force = new DenseVector(2) { [0] = 3, [1] = 3 }; + Vector force = new Vector(2) { [0] = 3, [1] = 3 }; NonLinearSolver Solver = NonLinearSolver.Builder .Solve(2, Reaction, Stiffness) .Under(force) - .WithInitialConditions(0.1, DenseVector.Create(2, 0), DenseVector.Create(2, 1)) + .WithInitialConditions(0.1, new Vector(2, 0), new Vector(2, 1)) .UsingArcLengthScheme(0.05) .NormalizeLoadWith(0.01) .Build(); @@ -184,9 +182,9 @@ static ILinearSolver Stiffness(Vector u) AssertSolutionsAreCorrect(Reaction, force, states); } - void AssertSolutionIsCorrect(Vector solution) + void AssertSolutionIsCorrect(Vector solution) { - double first = solution.First(); + double first = solution[0]; foreach (double d in solution) { Assert.Equal(first, d, decimalsPrecision); @@ -194,8 +192,8 @@ void AssertSolutionIsCorrect(Vector solution) } void AssertSolutionsAreCorrect( - Func, Vector> reaction - , Vector force + Func reaction + , Vector force , List states) { foreach (LoadState state in states) @@ -206,8 +204,8 @@ Func, Vector> reaction } void AssertAreCloseEnough( - Vector v1 - , Vector v2 + Vector v1 + , Vector v2 , double tolerance) => Assert.True((v1 - v2).Norm(2) <= tolerance); } diff --git a/NLES/Builder/SolverBuilder.cs b/NLES/Builder/SolverBuilder.cs index a059cff..d94af21 100644 --- a/NLES/Builder/SolverBuilder.cs +++ b/NLES/Builder/SolverBuilder.cs @@ -1,7 +1,5 @@ using System; -using MathNet.Numerics.LinearAlgebra; - namespace NLES { public partial class NonLinearSolver @@ -12,10 +10,10 @@ public class SolverBuilder public SolverBuilderStructure Solve( int degreesOfFreedom - , Func, Vector> reaction - , Func, ILinearSolver> stiffness) => new SolverBuilderStructure(degreesOfFreedom, Solver, reaction, stiffness); + , Func reaction + , Func stiffness) => new SolverBuilderStructure(degreesOfFreedom, Solver, reaction, stiffness); - public SolverBuilderIncrementalLoad Under(Vector referenceLoad) => new SolverBuilderIncrementalLoad(Solver, referenceLoad); + public SolverBuilderIncrementalLoad Under(Vector referenceLoad) => new SolverBuilderIncrementalLoad(Solver, referenceLoad); public SolverBuilderSchemeStandardNewtonRaphson UsingStandardNewtonRaphsonScheme(double loadFactorIncrement) => new SolverBuilderSchemeStandardNewtonRaphson(Solver, loadFactorIncrement); diff --git a/NLES/Builder/SolverBuilderIncrementalLoad.cs b/NLES/Builder/SolverBuilderIncrementalLoad.cs index e414ba0..138b318 100644 --- a/NLES/Builder/SolverBuilderIncrementalLoad.cs +++ b/NLES/Builder/SolverBuilderIncrementalLoad.cs @@ -1,5 +1,4 @@ -using MathNet.Numerics.LinearAlgebra; -using NLES.Contracts; +using NLES.Contracts; namespace NLES { @@ -7,13 +6,13 @@ public partial class NonLinearSolver { public class SolverBuilderIncrementalLoad : SolverBuilder { - public SolverBuilderIncrementalLoad(NonLinearSolver solver, Vector referenceLoad) + public SolverBuilderIncrementalLoad(NonLinearSolver solver, Vector referenceLoad) { Solver = solver; Solver.Info.ReferenceLoad = referenceLoad; } - public SolverBuilderIncrementalLoad WithInitialConditions(double lambda, Vector load, Vector displacement) + public SolverBuilderIncrementalLoad WithInitialConditions(double lambda, Vector load, Vector displacement) { Solver.Info.InitialLoad = load; Solver.State = new LoadState(lambda, displacement); diff --git a/NLES/Builder/SolverBuilderStructure.cs b/NLES/Builder/SolverBuilderStructure.cs index 22c9884..fb14760 100644 --- a/NLES/Builder/SolverBuilderStructure.cs +++ b/NLES/Builder/SolverBuilderStructure.cs @@ -1,6 +1,4 @@ -using MathNet.Numerics.LinearAlgebra; -using MathNet.Numerics.LinearAlgebra.Double; -using NLES.Contracts; +using NLES.Contracts; using System; namespace NLES @@ -12,26 +10,26 @@ public class SolverBuilderStructure : SolverBuilder public SolverBuilderStructure( int degreesOfFreedom , NonLinearSolver solver - , Func, Vector> reaction - , Func, ILinearSolver> stiffness) + , Func reaction + , Func stiffness) { CheckDregreesOfFreedom(degreesOfFreedom); CheckReaction(reaction); CheckStiffness(stiffness); Solver = solver; - Solver.State = new LoadState(0, new DenseVector(degreesOfFreedom)); + Solver.State = new LoadState(0, new Vector(degreesOfFreedom)); Solver.Info = new StructureInfo { - InitialLoad = new DenseVector(degreesOfFreedom), + InitialLoad = new Vector(degreesOfFreedom), Reaction = reaction, - ReferenceLoad = new DenseVector(degreesOfFreedom), + ReferenceLoad = new Vector(degreesOfFreedom), Stiffness = stiffness, }; Solver.Info.Reaction = reaction; } - private static void CheckStiffness(Func, ILinearSolver> stiffness) + private static void CheckStiffness(Func stiffness) { if (stiffness == null) { @@ -39,7 +37,7 @@ private static void CheckStiffness(Func, ILinearSolver> stiffness } } - private static void CheckReaction(Func, Vector> reaction) + private static void CheckReaction(Func reaction) { if (reaction == null) { diff --git a/NLES/Contracts/LoadState.cs b/NLES/Contracts/LoadState.cs index ef40b05..d0e51cd 100644 --- a/NLES/Contracts/LoadState.cs +++ b/NLES/Contracts/LoadState.cs @@ -1,13 +1,12 @@ -using MathNet.Numerics.LinearAlgebra; - + namespace NLES.Contracts { public struct LoadState { public double Lambda { get; set; } - public Vector Displacement { get; set; } + public Vector Displacement { get; set; } - internal LoadState(double lambda, Vector displacement) + internal LoadState(double lambda, Vector displacement) { Lambda = lambda; Displacement = displacement; diff --git a/NLES/Correction/CorrectionSchemeArcLength.cs b/NLES/Correction/CorrectionSchemeArcLength.cs index 3fc2065..7d3acc3 100644 --- a/NLES/Correction/CorrectionSchemeArcLength.cs +++ b/NLES/Correction/CorrectionSchemeArcLength.cs @@ -1,7 +1,6 @@ using System; using System.Collections.Generic; using System.Linq; -using MathNet.Numerics.LinearAlgebra; using NLES.Contracts; using NLES.Correction.Methods; @@ -18,8 +17,8 @@ internal class CorrectionSchemeArcLength : ICorrectionScheme public Result Correct(LoadState state, LoadIncrementalState prediction, StructureInfo info, - Vector dut, - Vector dur) + Vector dut, + Vector dur) { Result candidates = GetCandidates(dut, dur, prediction, info); @@ -34,8 +33,8 @@ public Result Correct(LoadState state, }; } - Result GetCandidates(Vector dut, - Vector dur, + Result GetCandidates(Vector dut, + Vector dur, LoadIncrementalState prediction, StructureInfo info) { diff --git a/NLES/Correction/CorrectionSchemeStandard.cs b/NLES/Correction/CorrectionSchemeStandard.cs index 0e725f4..407d235 100644 --- a/NLES/Correction/CorrectionSchemeStandard.cs +++ b/NLES/Correction/CorrectionSchemeStandard.cs @@ -1,6 +1,4 @@ -using MathNet.Numerics.LinearAlgebra; - -using NLES.Contracts; +using NLES.Contracts; namespace NLES.Correction { @@ -10,8 +8,8 @@ public Result Correct( LoadState state , LoadIncrementalState prediction , StructureInfo info - , Vector dut - , Vector dur) => new Result() + , Vector dut + , Vector dur) => new Result() { Value = new LoadIncrementalState(0, dur) }; diff --git a/NLES/Correction/Corrector.cs b/NLES/Correction/Corrector.cs index 530a0c4..e10fd42 100644 --- a/NLES/Correction/Corrector.cs +++ b/NLES/Correction/Corrector.cs @@ -2,8 +2,6 @@ using System.Collections.Generic; using System.Linq; -using MathNet.Numerics.LinearAlgebra; - using NLES.Contracts; namespace NLES.Correction @@ -56,9 +54,9 @@ LoadState state IEnumerable GetErrors(LoadState newState, LoadState oldState, StructureInfo info) { - Vector reaction = info.Reaction(newState.Displacement); - Vector equilibrium = info.InitialLoad + newState.Lambda * info.ReferenceLoad - reaction; - Vector incrementDisplacement = newState.Displacement - oldState.Displacement; + Vector reaction = info.Reaction(newState.Displacement); + Vector equilibrium = info.InitialLoad + newState.Lambda * info.ReferenceLoad - reaction; + Vector incrementDisplacement = newState.Displacement - oldState.Displacement; yield return incrementDisplacement.Norm(2) / newState.Displacement.Norm(2); yield return equilibrium.Norm(2) / info.ReferenceLoad.Norm(2); @@ -73,10 +71,10 @@ bool CheckConvergence(IEnumerable errors, IReadOnlyList toleranc Result GetCorrection(LoadState state, LoadIncrementalState prediction, StructureInfo info) { ILinearSolver stiffnessMatrix = info.Stiffness(state.Displacement); - Vector dut = stiffnessMatrix.Solve(info.ReferenceLoad); - Vector reaction = info.Reaction(state.Displacement); - Vector equilibrium = info.InitialLoad + state.Lambda * info.ReferenceLoad - reaction; - Vector dur = stiffnessMatrix.Solve(equilibrium); + Vector dut = stiffnessMatrix.Solve(info.ReferenceLoad); + Vector reaction = info.Reaction(state.Displacement); + Vector equilibrium = info.InitialLoad + state.Lambda * info.ReferenceLoad - reaction; + Vector dur = stiffnessMatrix.Solve(equilibrium); Result incrementalStateResult = Scheme.Correct(state, prediction, info, dut, dur); return incrementalStateResult; diff --git a/NLES/Correction/ICorrectionScheme.cs b/NLES/Correction/ICorrectionScheme.cs index 4a1b1e1..e48d002 100644 --- a/NLES/Correction/ICorrectionScheme.cs +++ b/NLES/Correction/ICorrectionScheme.cs @@ -1,5 +1,4 @@ -using MathNet.Numerics.LinearAlgebra; -using NLES.Contracts; +using NLES.Contracts; namespace NLES.Correction { @@ -8,7 +7,7 @@ internal interface ICorrectionScheme Result Correct(LoadState state, LoadIncrementalState prediction, StructureInfo info, - Vector dut, - Vector dur); + Vector dut, + Vector dur); } } \ No newline at end of file diff --git a/NLES/Correction/Methods/RestoringMethod.cs b/NLES/Correction/Methods/RestoringMethod.cs index 0ad7dac..49296be 100644 --- a/NLES/Correction/Methods/RestoringMethod.cs +++ b/NLES/Correction/Methods/RestoringMethod.cs @@ -1,5 +1,4 @@ -using MathNet.Numerics.LinearAlgebra; -using MoreLinq; +using MoreLinq; using NLES.Contracts; using System.Collections.Generic; using System.Linq; @@ -16,10 +15,10 @@ StructureInfo info { double Function(LoadIncrementalState candidate) { - Vector displacement = state.Displacement + candidate.IncrementDisplacement; - Vector reaction = info.Reaction(displacement); + Vector displacement = state.Displacement + candidate.IncrementDisplacement; + Vector reaction = info.Reaction(displacement); double lambda = state.Lambda + candidate.IncrementLambda; - Vector equilibriumVector = info.InitialLoad + lambda * info.ReferenceLoad - reaction; + Vector equilibriumVector = info.InitialLoad + lambda * info.ReferenceLoad - reaction; return equilibriumVector.Norm(2); } diff --git a/NLES/ILinearSolver.cs b/NLES/ILinearSolver.cs index 7c48fbe..f995cea 100644 --- a/NLES/ILinearSolver.cs +++ b/NLES/ILinearSolver.cs @@ -1,5 +1,4 @@ -using MathNet.Numerics.LinearAlgebra; - + namespace NLES { /// @@ -12,6 +11,6 @@ public interface ILinearSolver /// /// /// - Vector Solve(Vector input); + Vector Solve(Vector input); } } diff --git a/NLES/LoadIncrementalState.cs b/NLES/LoadIncrementalState.cs index 342d280..932810d 100644 --- a/NLES/LoadIncrementalState.cs +++ b/NLES/LoadIncrementalState.cs @@ -1,14 +1,13 @@ -using MathNet.Numerics.LinearAlgebra; - + namespace NLES { internal class LoadIncrementalState { internal double IncrementLambda { get; set; } - internal Vector IncrementDisplacement { get; set; } + internal Vector IncrementDisplacement { get; set; } - internal LoadIncrementalState(double incrementLambda, Vector incrementDisplacement) + internal LoadIncrementalState(double incrementLambda, Vector incrementDisplacement) { IncrementLambda = incrementLambda; IncrementDisplacement = incrementDisplacement; diff --git a/NLES/NLES.csproj b/NLES/NLES.csproj index f47f6ad..a865b19 100644 --- a/NLES/NLES.csproj +++ b/NLES/NLES.csproj @@ -5,7 +5,6 @@ - diff --git a/NLES/NonLinearSolver.cs b/NLES/NonLinearSolver.cs index 04a76bc..230cbd2 100644 --- a/NLES/NonLinearSolver.cs +++ b/NLES/NonLinearSolver.cs @@ -1,8 +1,6 @@ using System.Collections.Generic; using System.Runtime.CompilerServices; -using MathNet.Numerics.LinearAlgebra; - using NLES.Contracts; using NLES.Correction; using NLES.Prediction; @@ -22,7 +20,7 @@ public partial class NonLinearSolver public IEnumerable Broadcast() { ILinearSolver mK0 = Info.Stiffness(State.Displacement); - Vector Dv0 = mK0.Solve(Info.ReferenceLoad); + Vector Dv0 = mK0.Solve(Info.ReferenceLoad); double k0 = Info.ReferenceLoad.DotProduct(Dv0); while (true) { diff --git a/NLES/Prediction/IPredictionScheme.cs b/NLES/Prediction/IPredictionScheme.cs index cef159a..f54f95d 100644 --- a/NLES/Prediction/IPredictionScheme.cs +++ b/NLES/Prediction/IPredictionScheme.cs @@ -1,9 +1,8 @@ -using MathNet.Numerics.LinearAlgebra; - + namespace NLES.Prediction { internal interface IPredictionScheme { - double Predict(Vector Dvt, Vector fr); + double Predict(Vector Dvt, Vector fr); } } \ No newline at end of file diff --git a/NLES/Prediction/PredictionSchemeArcLength.cs b/NLES/Prediction/PredictionSchemeArcLength.cs index 071ea57..1b189a3 100644 --- a/NLES/Prediction/PredictionSchemeArcLength.cs +++ b/NLES/Prediction/PredictionSchemeArcLength.cs @@ -1,5 +1,4 @@ using System; -using MathNet.Numerics.LinearAlgebra; namespace NLES.Prediction { @@ -11,7 +10,7 @@ internal class PredictionSchemeArcLength : IPredictionScheme internal PredictionSchemeArcLength(double radius) => Radius = radius; - public double Predict(Vector Dvt, Vector fr) + public double Predict(Vector Dvt, Vector fr) { double dispDotProduct = Dvt.DotProduct(Dvt); double forceDotProduct = fr.DotProduct(fr); diff --git a/NLES/Prediction/PredictionSchemeStandard.cs b/NLES/Prediction/PredictionSchemeStandard.cs index 7a12d1b..d0750ff 100644 --- a/NLES/Prediction/PredictionSchemeStandard.cs +++ b/NLES/Prediction/PredictionSchemeStandard.cs @@ -1,5 +1,4 @@ -using MathNet.Numerics.LinearAlgebra; - + namespace NLES.Prediction { internal class PredictionSchemeStandard : IPredictionScheme @@ -8,6 +7,6 @@ internal class PredictionSchemeStandard : IPredictionScheme internal PredictionSchemeStandard(double dlambda) => DLambda = dlambda; - public double Predict(Vector Dvt, Vector fr) => DLambda; + public double Predict(Vector Dvt, Vector fr) => DLambda; } } diff --git a/NLES/Prediction/Predictor.cs b/NLES/Prediction/Predictor.cs index f380625..2664d99 100644 --- a/NLES/Prediction/Predictor.cs +++ b/NLES/Prediction/Predictor.cs @@ -1,5 +1,4 @@ -using MathNet.Numerics.LinearAlgebra; -using NLES.Contracts; +using NLES.Contracts; using System; namespace NLES.Prediction @@ -10,18 +9,18 @@ internal class Predictor internal LoadIncrementalState Predict(LoadState state, double initialStiffness, StructureInfo info) { - Vector equilibrium = info.InitialLoad + state.Lambda * info.ReferenceLoad - info.Reaction(state.Displacement); + Vector equilibrium = info.InitialLoad + state.Lambda * info.ReferenceLoad - info.Reaction(state.Displacement); ILinearSolver mK = info.Stiffness(state.Displacement); - Vector Dvt = mK.Solve(info.ReferenceLoad); - Vector Dvr = mK.Solve(equilibrium); + Vector Dvt = mK.Solve(info.ReferenceLoad); + Vector Dvr = mK.Solve(equilibrium); double bergam = GetBergamParameter(initialStiffness, Dvt, info); double DLambda = Scheme.Predict(Dvt, info.ReferenceLoad) * Math.Sign(bergam); - Vector Dv = DLambda * Dvt + Dvr; + Vector Dv = DLambda * Dvt + Dvr; return new LoadIncrementalState(DLambda, Dv); } - double GetBergamParameter(double k0, Vector Dvt, StructureInfo info) + double GetBergamParameter(double k0, Vector Dvt, StructureInfo info) => Math.Abs(k0 / info.ReferenceLoad.DotProduct(Dvt)); } } diff --git a/NLES/StructureInfo.cs b/NLES/StructureInfo.cs index 4788d71..8b50e56 100644 --- a/NLES/StructureInfo.cs +++ b/NLES/StructureInfo.cs @@ -1,14 +1,12 @@ using System; -using MathNet.Numerics.LinearAlgebra; -using MathNet.Numerics.LinearAlgebra.Double; namespace NLES { internal class StructureInfo { - internal Func, Vector> Reaction { get; set; } - internal Vector ReferenceLoad { get; set; } - internal Vector InitialLoad { get; set; } - internal Func, ILinearSolver> Stiffness { get; set; } + internal Func Reaction { get; set; } + internal Vector ReferenceLoad { get; set; } + internal Vector InitialLoad { get; set; } + internal Func Stiffness { get; set; } } } diff --git a/NLES/Vector.cs b/NLES/Vector.cs new file mode 100644 index 0000000..334c855 --- /dev/null +++ b/NLES/Vector.cs @@ -0,0 +1,88 @@ +using System; +using System.Collections; + +namespace NLES +{ + public class Vector : IEnumerable + { + private readonly double[] array; + + public Vector(int dimension, double value = 0) + { + array = new double[dimension]; + for (int i = 0; i < dimension; i++) + { + array[i] = value; + } + } + + public int Dimension => array.Length; + + public double this[int index] + { + get => array[index]; + set => array[index] = value; + } + + public double DotProduct(Vector vector) + { + double result = 0; + for (int i = 0; i < Dimension; i++) + { + result += this[i] * vector[i]; + } + + return result; + } + + public double Norm(double exp) + { + double result = 0; + for (int i = 0; i < Dimension; i++) + { + result += Math.Pow(this[i], exp); + } + result = Math.Pow(result, 1.0 / exp); + + return result; + } + + public IEnumerator GetEnumerator() => array.GetEnumerator(); + + public static Vector operator +(Vector v1, Vector v2) + { + Vector result = new Vector(v1.Dimension); + + for (int i = 0; i < v1.Dimension; i++) + { + result[i] = v1[i] + v2[i]; + } + + return result; + } + + public static Vector operator -(Vector v1, Vector v2) + { + Vector result = new Vector(v1.Dimension); + + for (int i = 0; i < v1.Dimension; i++) + { + result[i] = v1[i] - v2[i]; + } + + return result; + } + + public static Vector operator *(double factor, Vector v) + { + Vector result = new Vector(v.Dimension); + + for (int i = 0; i < v.Dimension; i++) + { + result[i] = factor * v[i]; + } + + return result; + } + } +}