From 9b1529f3ed16660a2364edff6b656dfa90ce1b46 Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Thu, 30 Nov 2023 20:58:49 -0700 Subject: [PATCH 01/12] Templates working again --- OndselSolver/ASMTAssembly.cpp | 20 +- OndselSolver/ASMTPrincipalMassMarker.cpp | 6 +- OndselSolver/ASMTPrincipalMassMarker.h | 2 +- OndselSolver/ASMTSpatialItem.cpp | 4 +- OndselSolver/ASMTSpatialItem.h | 2 +- OndselSolver/AngleZIeqcJec.cpp | 4 +- OndselSolver/AngleZIeqcJeqc.cpp | 8 +- OndselSolver/AnyPosICNewtonRaphson.cpp | 2 +- OndselSolver/CADSystem.cpp | 86 +-- OndselSolver/ConstVelConstraintIqcJc.cpp | 4 +- OndselSolver/ConstVelConstraintIqcJqc.cpp | 8 +- OndselSolver/DiagonalMatrix.cpp | 87 +-- OndselSolver/DiagonalMatrix.h | 133 +++- OndselSolver/DiagonalMatrix.ref.h | 9 - OndselSolver/DifferenceOperator.cpp | 2 +- OndselSolver/DirectionCosineIecJec.cpp | 2 +- OndselSolver/DirectionCosineIeqcJec.cpp | 4 +- OndselSolver/DirectionCosineIeqcJeqc.cpp | 8 +- OndselSolver/DirectionCosineIeqctJeqc.cpp | 6 +- OndselSolver/DispCompIecJecIe.cpp | 2 +- OndselSolver/DispCompIecJecKeqc.cpp | 8 +- OndselSolver/DispCompIeqcJecIe.cpp | 14 +- OndselSolver/DispCompIeqcJecKeqc.cpp | 12 +- OndselSolver/DispCompIeqcJeqcIe.cpp | 16 +- OndselSolver/DispCompIeqcJeqcKeqc.cpp | 12 +- OndselSolver/DispCompIeqcJeqcKeqct.cpp | 8 +- OndselSolver/DispCompIeqctJeqcIe.cpp | 12 +- OndselSolver/DispCompIeqctJeqcKeqct.cpp | 6 +- OndselSolver/DistIeqcJec.cpp | 12 +- OndselSolver/DistIeqcJeqc.cpp | 16 +- OndselSolver/DistxyIeqcJec.cpp | 6 +- OndselSolver/DistxyIeqcJeqc.cpp | 14 +- OndselSolver/EndFramec.h | 2 +- OndselSolver/EndFrameqc.cpp | 12 +- OndselSolver/EndFrameqct.cpp | 655 +++++++++-------- OndselSolver/EndFrameqct2.cpp | 3 +- OndselSolver/EulerAngles.cpp | 73 -- OndselSolver/EulerAngles.h | 90 ++- OndselSolver/EulerAnglesDDot.h | 6 +- OndselSolver/EulerAnglesDot.h | 6 +- OndselSolver/EulerAngleszxz.h | 6 +- OndselSolver/EulerAngleszxzDDot.h | 12 +- OndselSolver/EulerAngleszxzDot.h | 6 +- OndselSolver/EulerParameters.cpp | 305 +------- OndselSolver/EulerParameters.h | 303 +++++++- OndselSolver/EulerParametersDot.h | 8 +- OndselSolver/FullColumn.cpp | 191 +---- OndselSolver/FullColumn.h | 239 +++++- OndselSolver/FullColumn.ref.h | 14 - OndselSolver/FullMatrix.cpp | 674 +---------------- OndselSolver/FullMatrix.h | 844 ++++++++++++++++++---- OndselSolver/FullMatrix.ref.h | 19 - OndselSolver/FullRow.cpp | 51 +- OndselSolver/FullRow.h | 70 +- OndselSolver/FullRow.ref.h | 15 - OndselSolver/GeneralSpline.cpp | 2 +- OndselSolver/Item.h | 4 - OndselSolver/LDUFullMat.cpp | 2 +- OndselSolver/MBDynBody.cpp | 6 +- OndselSolver/MBDynItem.cpp | 10 +- OndselSolver/MBDynJoint.cpp | 2 +- OndselSolver/MBDynMarker.cpp | 4 +- OndselSolver/MBDynStructural.cpp | 2 +- OndselSolver/MarkerFrame.cpp | 2 +- OndselSolver/MarkerFrame.h | 4 +- OndselSolver/MomentOfInertiaSolver.cpp | 26 +- OndselSolver/OndselSolver.vcxproj | 4 - OndselSolver/OndselSolver.vcxproj.filters | 12 - OndselSolver/OrbitAngleZIeqcJec.cpp | 6 +- OndselSolver/OrbitAngleZIeqcJeqc.cpp | 14 +- OndselSolver/Part.cpp | 16 +- OndselSolver/Part.h | 3 - OndselSolver/RackPinConstraintIqcJc.cpp | 4 +- OndselSolver/RackPinConstraintIqcJqc.cpp | 6 +- OndselSolver/RowTypeMatrix.h | 1 - OndselSolver/ScrewConstraintIqcJc.cpp | 4 +- OndselSolver/ScrewConstraintIqcJqc.cpp | 6 +- OndselSolver/StableBackwardDifference.cpp | 4 +- OndselSolver/VelICSolver.cpp | 2 +- 79 files changed, 2063 insertions(+), 2222 deletions(-) delete mode 100644 OndselSolver/DiagonalMatrix.ref.h delete mode 100644 OndselSolver/FullColumn.ref.h delete mode 100644 OndselSolver/FullMatrix.ref.h delete mode 100644 OndselSolver/FullRow.ref.h diff --git a/OndselSolver/ASMTAssembly.cpp b/OndselSolver/ASMTAssembly.cpp index c0929d2..0897438 100644 --- a/OndselSolver/ASMTAssembly.cpp +++ b/OndselSolver/ASMTAssembly.cpp @@ -263,7 +263,7 @@ void MbD::ASMTAssembly::runSinglePendulum() assembly->setName(str); auto pos3D = std::make_shared>(ListD{ 0, 0, 0 }); assembly->setPosition3D(pos3D); - auto rotMat = std::make_shared(ListListD{ + auto rotMat = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -277,11 +277,11 @@ void MbD::ASMTAssembly::runSinglePendulum() auto massMarker = std::make_shared(); massMarker->setMass(0.0); massMarker->setDensity(0.0); - auto aJ = std::make_shared(ListD{ 0, 0, 0 }); + auto aJ = std::make_shared>(ListD{ 0, 0, 0 }); massMarker->setMomentOfInertias(aJ); pos3D = std::make_shared>(ListD{ 0, 0, 0 }); massMarker->setPosition3D(pos3D); - rotMat = std::make_shared(ListListD{ + rotMat = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -294,7 +294,7 @@ void MbD::ASMTAssembly::runSinglePendulum() mkr->setName(str); pos3D = std::make_shared>(ListD{ 0, 0, 0 }); mkr->setPosition3D(pos3D); - rotMat = std::make_shared(ListListD{ + rotMat = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -307,7 +307,7 @@ void MbD::ASMTAssembly::runSinglePendulum() part->setName(str); pos3D = std::make_shared>(ListD{ -0.1, -0.1, -0.1 }); part->setPosition3D(pos3D); - rotMat = std::make_shared(ListListD{ + rotMat = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -322,11 +322,11 @@ void MbD::ASMTAssembly::runSinglePendulum() massMarker = std::make_shared(); massMarker->setMass(0.2); massMarker->setDensity(10.0); - aJ = std::make_shared(ListD{ 8.3333333333333e-4, 0.016833333333333, 0.017333333333333 }); + aJ = std::make_shared>(ListD{ 8.3333333333333e-4, 0.016833333333333, 0.017333333333333 }); massMarker->setMomentOfInertias(aJ); pos3D = std::make_shared>(ListD{ 0.5, 0.1, 0.05 }); massMarker->setPosition3D(pos3D); - rotMat = std::make_shared(ListListD{ + rotMat = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -339,7 +339,7 @@ void MbD::ASMTAssembly::runSinglePendulum() mkr->setName(str); pos3D = std::make_shared>(ListD{ 0.1, 0.1, 0.1 }); mkr->setPosition3D(pos3D); - rotMat = std::make_shared(ListListD{ + rotMat = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -1052,9 +1052,9 @@ void MbD::ASMTAssembly::initprincipalMassMarker() principalMassMarker = std::make_shared(); principalMassMarker->mass = 0.0; principalMassMarker->density = 0.0; - principalMassMarker->momentOfInertias = std::make_shared(3, 0); + principalMassMarker->momentOfInertias = std::make_shared>(3, 0); //principalMassMarker->position3D = std::make_shared>(3, 0); - //principalMassMarker->rotationMatrix = FullMatrixDouble>::identitysptr(3); + //principalMassMarker->rotationMatrix = FullMatrix>::identitysptr(3); } std::shared_ptr MbD::ASMTAssembly::spatialContainerAt(std::shared_ptr self, std::string& longname) diff --git a/OndselSolver/ASMTPrincipalMassMarker.cpp b/OndselSolver/ASMTPrincipalMassMarker.cpp index 1e4e838..003c4de 100644 --- a/OndselSolver/ASMTPrincipalMassMarker.cpp +++ b/OndselSolver/ASMTPrincipalMassMarker.cpp @@ -31,7 +31,7 @@ void MbD::ASMTPrincipalMassMarker::parseASMT(std::vector& lines) lines.erase(lines.begin()); assert(lines[0] == (leadingTabs + "RotationMatrix")); lines.erase(lines.begin()); - rotationMatrix = std::make_shared(3); + rotationMatrix = std::make_shared>(3); for (int i = 0; i < 3; i++) { auto row = readRowOfDoubles(lines[0]); @@ -44,7 +44,7 @@ void MbD::ASMTPrincipalMassMarker::parseASMT(std::vector& lines) lines.erase(lines.begin()); assert(lines[0] == (leadingTabs + "MomentOfInertias")); lines.erase(lines.begin()); - momentOfInertias = std::make_shared(3); + momentOfInertias = std::make_shared>(3); auto row = readRowOfDoubles(lines[0]); lines.erase(lines.begin()); for (int i = 0; i < 3; i++) @@ -75,7 +75,7 @@ void MbD::ASMTPrincipalMassMarker::setMomentOfInertias(DiagMatDsptr mat) // Overloads to simplify syntax. void MbD::ASMTPrincipalMassMarker::setMomentOfInertias(double a, double b, double c) { - momentOfInertias = std::make_shared(ListD{ a, b, c }); + momentOfInertias = std::make_shared>(ListD{ a, b, c }); } void MbD::ASMTPrincipalMassMarker::storeOnLevel(std::ofstream& os, int level) diff --git a/OndselSolver/ASMTPrincipalMassMarker.h b/OndselSolver/ASMTPrincipalMassMarker.h index 633ad52..c4eed88 100644 --- a/OndselSolver/ASMTPrincipalMassMarker.h +++ b/OndselSolver/ASMTPrincipalMassMarker.h @@ -27,7 +27,7 @@ namespace MbD { double mass = 1.0; double density = 10.0; - DiagMatDsptr momentOfInertias = std::make_shared(ListD{ 1.0, 2.0, 3.0 }); + DiagMatDsptr momentOfInertias = std::make_shared>(ListD{ 1.0, 2.0, 3.0 }); }; } diff --git a/OndselSolver/ASMTSpatialItem.cpp b/OndselSolver/ASMTSpatialItem.cpp index 21c920f..b64f8f8 100644 --- a/OndselSolver/ASMTSpatialItem.cpp +++ b/OndselSolver/ASMTSpatialItem.cpp @@ -48,7 +48,7 @@ void MbD::ASMTSpatialItem::readRotationMatrix(std::vector& lines) { assert(lines[0].find("RotationMatrix") != std::string::npos); lines.erase(lines.begin()); - rotationMatrix = std::make_shared(3, 0); + rotationMatrix = std::make_shared>(3, 0); for (int i = 0; i < 3; i++) { auto& row = rotationMatrix->at(i); @@ -87,7 +87,7 @@ void MbD::ASMTSpatialItem::setRotationMatrix(double v11, double v12, double v13, double v21, double v22, double v23, double v31, double v32, double v33) { - rotationMatrix = std::make_shared(ListListD{ + rotationMatrix = std::make_shared>(ListListD{ { v11, v12, v13 }, { v21, v22, v23 }, { v31, v32, v33 } diff --git a/OndselSolver/ASMTSpatialItem.h b/OndselSolver/ASMTSpatialItem.h index 0469336..29bb7ad 100644 --- a/OndselSolver/ASMTSpatialItem.h +++ b/OndselSolver/ASMTSpatialItem.h @@ -35,7 +35,7 @@ namespace MbD { FMatDsptr getRotationMatrix(int i); FColDsptr position3D = std::make_shared>(3); - FMatDsptr rotationMatrix = std::make_shared(ListListD{ + FMatDsptr rotationMatrix = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } diff --git a/OndselSolver/AngleZIeqcJec.cpp b/OndselSolver/AngleZIeqcJec.cpp index e5b4de6..e315ea3 100644 --- a/OndselSolver/AngleZIeqcJec.cpp +++ b/OndselSolver/AngleZIeqcJec.cpp @@ -18,7 +18,7 @@ MbD::AngleZIeqcJec::AngleZIeqcJec() MbD::AngleZIeqcJec::AngleZIeqcJec(EndFrmsptr frmi, EndFrmsptr frmj) : AngleZIecJec(frmi, frmj) { pthezpEI = std::make_shared>(4); - ppthezpEIpEI = std::make_shared(4, 4); + ppthezpEIpEI = std::make_shared>(4, 4); } void MbD::AngleZIeqcJec::calcPostDynCorrectorIteration() @@ -67,7 +67,7 @@ void MbD::AngleZIeqcJec::initialize() { AngleZIecJec::initialize(); pthezpEI = std::make_shared>(4); - ppthezpEIpEI = std::make_shared(4, 4); + ppthezpEIpEI = std::make_shared>(4, 4); } FMatDsptr MbD::AngleZIeqcJec::ppvaluepEIpEI() diff --git a/OndselSolver/AngleZIeqcJeqc.cpp b/OndselSolver/AngleZIeqcJeqc.cpp index 0ca4694..6b50183 100644 --- a/OndselSolver/AngleZIeqcJeqc.cpp +++ b/OndselSolver/AngleZIeqcJeqc.cpp @@ -19,8 +19,8 @@ MbD::AngleZIeqcJeqc::AngleZIeqcJeqc() MbD::AngleZIeqcJeqc::AngleZIeqcJeqc(EndFrmsptr frmi, EndFrmsptr frmj) : AngleZIeqcJec(frmi, frmj) { pthezpEJ = std::make_shared>(4); - ppthezpEIpEJ = std::make_shared(4, 4); - ppthezpEJpEJ = std::make_shared(4, 4); + ppthezpEIpEJ = std::make_shared>(4, 4); + ppthezpEJpEJ = std::make_shared>(4, 4); } void MbD::AngleZIeqcJeqc::calcPostDynCorrectorIteration() @@ -88,8 +88,8 @@ void MbD::AngleZIeqcJeqc::initialize() { AngleZIeqcJec::initialize(); pthezpEJ = std::make_shared>(4); - ppthezpEIpEJ = std::make_shared(4, 4); - ppthezpEJpEJ = std::make_shared(4, 4); + ppthezpEIpEJ = std::make_shared>(4, 4); + ppthezpEJpEJ = std::make_shared>(4, 4); } FMatDsptr MbD::AngleZIeqcJeqc::ppvaluepEIpEJ() diff --git a/OndselSolver/AnyPosICNewtonRaphson.cpp b/OndselSolver/AnyPosICNewtonRaphson.cpp index 7e4fd46..824a0e1 100644 --- a/OndselSolver/AnyPosICNewtonRaphson.cpp +++ b/OndselSolver/AnyPosICNewtonRaphson.cpp @@ -32,7 +32,7 @@ void AnyPosICNewtonRaphson::initializeGlobally() void AnyPosICNewtonRaphson::createVectorsAndMatrices() { qsuOld = std::make_shared>(nqsu); - qsuWeights = std::make_shared(nqsu); + qsuWeights = std::make_shared>(nqsu); SystemNewtonRaphson::createVectorsAndMatrices(); } diff --git a/OndselSolver/CADSystem.cpp b/OndselSolver/CADSystem.cpp index 9737bba..a3e00ff 100644 --- a/OndselSolver/CADSystem.cpp +++ b/OndselSolver/CADSystem.cpp @@ -87,9 +87,9 @@ void CADSystem::runOndselSinglePendulum() auto assembly1 = CREATE::With("/Assembly1"); std::cout << "assembly1->name " << assembly1->name << std::endl; assembly1->m = 0.0; - assembly1->aJ = std::make_shared(ListD{ 0, 0, 0 }); + assembly1->aJ = std::make_shared>(ListD{ 0, 0, 0 }); qX = std::make_shared>(ListD{ 0, 0, 0 }); - aAap = std::make_shared(ListListD{ + aAap = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -106,7 +106,7 @@ void CADSystem::runOndselSinglePendulum() auto marker2 = CREATE::With("/Assembly1/Marker2"); rpmp = std::make_shared>(ListD{ 0.0, 0.0, 0.0 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -117,7 +117,7 @@ void CADSystem::runOndselSinglePendulum() auto marker1 = CREATE::With("/Assembly1/Marker1"); rpmp = std::make_shared>(ListD{ 0.0, 3.0, 0.0 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 0, 1 }, { 0, -1, 0 } @@ -130,9 +130,9 @@ void CADSystem::runOndselSinglePendulum() auto crankPart1 = CREATE::With("/Assembly1/Part1"); std::cout << "crankPart1->name " << crankPart1->name << std::endl; crankPart1->m = 1.0; - crankPart1->aJ = std::make_shared(ListD{ 1, 1, 1 }); + crankPart1->aJ = std::make_shared>(ListD{ 1, 1, 1 }); qX = std::make_shared>(ListD{ 0.4, 0.0, -0.05 }); - aAap = std::make_shared(ListListD{ + aAap = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -149,7 +149,7 @@ void CADSystem::runOndselSinglePendulum() auto marker1 = CREATE::With("/Assembly1/Part1/Marker1"); rpmp = std::make_shared>(ListD{ -0.4, 0.0, 0.05 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -160,7 +160,7 @@ void CADSystem::runOndselSinglePendulum() auto marker2 = CREATE::With("/Assembly1/Part1/Marker2"); rpmp = std::make_shared>(ListD{ 0.4, 0.0, 0.05 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -222,9 +222,9 @@ void CADSystem::runOndselDoublePendulum() auto assembly1 = CREATE::With("/Assembly1"); std::cout << "assembly1->name " << assembly1->name << std::endl; assembly1->m = 0.0; - assembly1->aJ = std::make_shared(ListD{ 0, 0, 0 }); + assembly1->aJ = std::make_shared>(ListD{ 0, 0, 0 }); qX = std::make_shared>(ListD{ 0, 0, 0 }); - aAap = std::make_shared(ListListD{ + aAap = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -241,7 +241,7 @@ void CADSystem::runOndselDoublePendulum() auto marker2 = CREATE::With("/Assembly1/Marker2"); rpmp = std::make_shared>(ListD{ 0.0, 0.0, 0.0 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -252,7 +252,7 @@ void CADSystem::runOndselDoublePendulum() auto marker1 = CREATE::With("/Assembly1/Marker1"); rpmp = std::make_shared>(ListD{ 0.0, 3.0, 0.0 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 0, 1 }, { 0, -1, 0 } @@ -265,9 +265,9 @@ void CADSystem::runOndselDoublePendulum() auto crankPart1 = CREATE::With("/Assembly1/Part1"); std::cout << "crankPart1->name " << crankPart1->name << std::endl; crankPart1->m = 1.0; - crankPart1->aJ = std::make_shared(ListD{ 1, 1, 1 }); + crankPart1->aJ = std::make_shared>(ListD{ 1, 1, 1 }); qX = std::make_shared>(ListD{ 0.4, 0.0, -0.05 }); - aAap = std::make_shared(ListListD{ + aAap = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -284,7 +284,7 @@ void CADSystem::runOndselDoublePendulum() auto marker1 = CREATE::With("/Assembly1/Part1/Marker1"); rpmp = std::make_shared>(ListD{ -0.4, 0.0, 0.05 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -295,7 +295,7 @@ void CADSystem::runOndselDoublePendulum() auto marker2 = CREATE::With("/Assembly1/Part1/Marker2"); rpmp = std::make_shared>(ListD{ 0.4, 0.0, 0.05 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -307,7 +307,7 @@ void CADSystem::runOndselDoublePendulum() auto conrodPart2 = CREATE::With("/Assembly1/Part2"); std::cout << "conrodPart2->name " << conrodPart2->name << std::endl; conrodPart2->m = 1.0; - conrodPart2->aJ = std::make_shared(ListD{ 1, 1, 1 }); + conrodPart2->aJ = std::make_shared>(ListD{ 1, 1, 1 }); qX = std::make_shared>(ListD{ 0.15, 0.1, 0.05 }); qE = std::make_shared>(ListD{ 0.0, 0.0, 1.0, 0.0 }); auto eulerParameters = CREATE>::With(ListD{ 0.0, 0.0, 1.0, 0.0 }); @@ -325,7 +325,7 @@ void CADSystem::runOndselDoublePendulum() auto marker1 = CREATE::With("/Assembly1/Part2/Marker1"); rpmp = std::make_shared>(ListD{ -0.65, 0.0, -0.05 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0} @@ -336,7 +336,7 @@ void CADSystem::runOndselDoublePendulum() auto marker2 = CREATE::With("/Assembly1/Part2/Marker2"); rpmp = std::make_shared>(ListD{ 0.65, 0.0, -0.05 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0} @@ -396,7 +396,7 @@ void CADSystem::runOndselPiston() auto assembly1 = CREATE::With("/Assembly1"); std::cout << "assembly1->name " << assembly1->name << std::endl; assembly1->m = 0.0; - assembly1->aJ = std::make_shared(ListD{ 0, 0, 0 }); + assembly1->aJ = std::make_shared>(ListD{ 0, 0, 0 }); qX = std::make_shared>(ListD{ 0, 0, 0 }); qE = std::make_shared>(ListD{ 0, 0, 0, 1 }); assembly1->setqX(qX); @@ -417,7 +417,7 @@ void CADSystem::runOndselPiston() auto marker2 = CREATE::With("/Assembly1/Marker2"); rpmp = std::make_shared>(ListD{ 0.0, 0.0, 0.0 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -428,7 +428,7 @@ void CADSystem::runOndselPiston() auto marker1 = CREATE::With("/Assembly1/Marker1"); rpmp = std::make_shared>(ListD{ 0.0, 3.0, 0.0 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 0, 1 }, { 0, -1, 0 } @@ -441,7 +441,7 @@ void CADSystem::runOndselPiston() auto crankPart1 = CREATE::With("/Assembly1/Part1"); std::cout << "crankPart1->name " << crankPart1->name << std::endl; crankPart1->m = 1.0; - crankPart1->aJ = std::make_shared(ListD{ 1, 1, 1 }); + crankPart1->aJ = std::make_shared>(ListD{ 1, 1, 1 }); qX = std::make_shared>(ListD{ 0.4, 0.0, -0.05 }); qE = std::make_shared>(ListD{ 0.0, 0.0, 0.0, 1.0 }); crankPart1->setqX(qX); @@ -460,7 +460,7 @@ void CADSystem::runOndselPiston() auto marker1 = CREATE::With("/Assembly1/Part1/Marker1"); rpmp = std::make_shared>(ListD{ -0.4, 0.0, 0.05 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -471,7 +471,7 @@ void CADSystem::runOndselPiston() auto marker2 = CREATE::With("/Assembly1/Part1/Marker2"); rpmp = std::make_shared>(ListD{ 0.4, 0.0, 0.05 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -483,7 +483,7 @@ void CADSystem::runOndselPiston() auto conrodPart2 = CREATE::With("/Assembly1/Part2"); std::cout << "conrodPart2->name " << conrodPart2->name << std::endl; conrodPart2->m = 1.0; - conrodPart2->aJ = std::make_shared(ListD{ 1, 1, 1 }); + conrodPart2->aJ = std::make_shared>(ListD{ 1, 1, 1 }); qX = std::make_shared>(ListD{ 0.15, 0.1, 0.05 }); qE = std::make_shared>(ListD{ 0.0, 0.0, 1.0, 0.0 }); conrodPart2->setqX(qX); @@ -502,7 +502,7 @@ void CADSystem::runOndselPiston() auto marker1 = CREATE::With("/Assembly1/Part2/Marker1"); rpmp = std::make_shared>(ListD{ -0.65, 0.0, -0.05 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0} @@ -513,7 +513,7 @@ void CADSystem::runOndselPiston() auto marker2 = CREATE::With("/Assembly1/Part2/Marker2"); rpmp = std::make_shared>(ListD{ 0.65, 0.0, -0.05 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0} @@ -525,7 +525,7 @@ void CADSystem::runOndselPiston() auto pistonPart3 = CREATE::With("/Assembly1/Part3"); std::cout << "pistonPart3->name " << pistonPart3->name << std::endl; pistonPart3->m = 1.0; - pistonPart3->aJ = std::make_shared(ListD{ 1, 1, 1 }); + pistonPart3->aJ = std::make_shared>(ListD{ 1, 1, 1 }); qX = std::make_shared>(ListD{ -0.0, 1.5, 0.0 }); qE = std::make_shared>(ListD{ 0.70710678118655, 0.70710678118655, 0.0, 0.0 }); pistonPart3->setqX(qX); @@ -544,7 +544,7 @@ void CADSystem::runOndselPiston() auto marker1 = CREATE::With("/Assembly1/Part3/Marker1"); rpmp = std::make_shared>(ListD{ -0.5, 0.0, 0.0 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ {0.0, 1.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 0.0, -1.0} @@ -555,7 +555,7 @@ void CADSystem::runOndselPiston() auto marker2 = CREATE::With("/Assembly1/Part3/Marker2"); rpmp = std::make_shared>(ListD{ 0.5, 0.0, 0.0 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0} @@ -633,7 +633,7 @@ void CADSystem::runPiston() auto assembly1 = CREATE::With("/Assembly1"); std::cout << "assembly1->name " << assembly1->name << std::endl; assembly1->m = 0.0; - assembly1->aJ = std::make_shared(ListD{ 0, 0, 0 }); + assembly1->aJ = std::make_shared>(ListD{ 0, 0, 0 }); qX = std::make_shared>(ListD{ 0, 0, 0 }); qE = std::make_shared>(ListD{ 0, 0, 0, 1 }); assembly1->setqX(qX); @@ -654,7 +654,7 @@ void CADSystem::runPiston() auto marker2 = CREATE::With("/Assembly1/Marker2"); rpmp = std::make_shared>(ListD{ 0.0, 0.0, 0.0 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -665,7 +665,7 @@ void CADSystem::runPiston() auto marker1 = CREATE::With("/Assembly1/Marker1"); rpmp = std::make_shared>(ListD{ 0.0, 2.8817526385684, 0.0 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 0, 1 }, { 0, -1, 0 } @@ -678,7 +678,7 @@ void CADSystem::runPiston() auto crankPart1 = CREATE::With("/Assembly1/Part1"); std::cout << "crankPart1->name " << crankPart1->name << std::endl; crankPart1->m = 0.045210530089461; - crankPart1->aJ = std::make_shared(ListD{ 1.7381980042084e-4, 0.003511159968501, 0.0036154518487535 }); + crankPart1->aJ = std::make_shared>(ListD{ 1.7381980042084e-4, 0.003511159968501, 0.0036154518487535 }); qX = std::make_shared>(ListD{ 0.38423368514246, 2.6661567755108e-17, -0.048029210642807 }); qE = std::make_shared>(ListD{ 0.0, 0.0, 0.0, 1.0 }); crankPart1->setqX(qX); @@ -697,7 +697,7 @@ void CADSystem::runPiston() auto marker1 = CREATE::With("/Assembly1/Part1/Marker1"); rpmp = std::make_shared>(ListD{ -0.38423368514246, -2.6661567755108e-17, 0.048029210642807 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -708,7 +708,7 @@ void CADSystem::runPiston() auto marker2 = CREATE::With("/Assembly1/Part1/Marker2"); rpmp = std::make_shared>(ListD{ 0.38423368514246, -2.6661567755108e-17, 0.048029210642807 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } @@ -720,7 +720,7 @@ void CADSystem::runPiston() auto conrodPart2 = CREATE::With("/Assembly1/Part2"); std::cout << "conrodPart2->name " << conrodPart2->name << std::endl; conrodPart2->m = 0.067815795134192; - conrodPart2->aJ = std::make_shared(ListD{ 2.6072970063126e-4, 0.011784982468533, 0.011941420288912 }); + conrodPart2->aJ = std::make_shared>(ListD{ 2.6072970063126e-4, 0.011784982468533, 0.011941420288912 }); qX = std::make_shared>(ListD{ 0.38423368514246, 0.49215295678475, 0.048029210642807 }); qE = std::make_shared>(ListD{ 0.0, 0.0, 0.89871703427292, 0.43852900965351 }); conrodPart2->setqX(qX); @@ -739,7 +739,7 @@ void CADSystem::runPiston() auto marker1 = CREATE::With("/Assembly1/Part2/Marker1"); rpmp = std::make_shared>(ListD{ -0.6243797383565, 1.1997705489799e-16, -0.048029210642807 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ {1.0, 2.7755575615629e-16, 0.0}, {-2.7755575615629e-16, 1.0, 0.0}, {0.0, 0.0, 1.0} @@ -750,7 +750,7 @@ void CADSystem::runPiston() auto marker2 = CREATE::With("/Assembly1/Part2/Marker2"); rpmp = std::make_shared>(ListD{ 0.6243797383565, -2.1329254204087e-16, -0.048029210642807 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ {1.0, 2.4980018054066e-16, 2.2204460492503e-16}, {-2.4980018054066e-16, 1.0, 4.1633363423443e-17}, {-2.2204460492503e-16, -4.1633363423443e-17, 1.0} @@ -762,7 +762,7 @@ void CADSystem::runPiston() auto pistonPart3 = CREATE::With("/Assembly1/Part3"); std::cout << "pistonPart3->name " << pistonPart3->name << std::endl; pistonPart3->m = 1.730132083368; - pistonPart3->aJ = std::make_shared(ListD{ 0.19449049546716, 0.23028116340971, 0.23028116340971 }); + pistonPart3->aJ = std::make_shared>(ListD{ 0.19449049546716, 0.23028116340971, 0.23028116340971 }); qX = std::make_shared>(ListD{ -1.283972762056e-18, 1.4645980199976, -4.7652385308244e-17 }); qE = std::make_shared>(ListD{ 0.70710678118655, 0.70710678118655, 0.0, 0.0 }); pistonPart3->setqX(qX); @@ -781,7 +781,7 @@ void CADSystem::runPiston() auto marker1 = CREATE::With("/Assembly1/Part3/Marker1"); rpmp = std::make_shared>(ListD{ -0.48029210642807, 7.6201599718927e-18, -2.816737703896e-17 }); marker1->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ {9.2444637330587e-33, 1.0, 2.2204460492503e-16}, {1.0, -9.2444637330587e-33, -1.0785207688569e-32}, {-1.0785207688569e-32, 2.2204460492503e-16, -1.0} @@ -792,7 +792,7 @@ void CADSystem::runPiston() auto marker2 = CREATE::With("/Assembly1/Part3/Marker2"); rpmp = std::make_shared>(ListD{ 0.48029210642807, 1.7618247880058e-17, 2.5155758471256e-17 }); marker2->setrpmp(rpmp); - aApm = std::make_shared(ListListD{ + aApm = std::make_shared>(ListListD{ {6.9388939039072e-18, -6.4146353042213e-50, 1.0}, {1.0, -6.9388939039072e-18, 6.9388939039072e-18}, {-6.9388939039072e-18, 1.0, -7.4837411882581e-50} diff --git a/OndselSolver/ConstVelConstraintIqcJc.cpp b/OndselSolver/ConstVelConstraintIqcJc.cpp index 117c2c2..22379fe 100644 --- a/OndselSolver/ConstVelConstraintIqcJc.cpp +++ b/OndselSolver/ConstVelConstraintIqcJc.cpp @@ -19,7 +19,7 @@ using namespace MbD; MbD::ConstVelConstraintIqcJc::ConstVelConstraintIqcJc(EndFrmsptr frmi, EndFrmsptr frmj) : ConstVelConstraintIJ(frmi, frmj) { pGpEI = std::make_shared>(4); - ppGpEIpEI = std::make_shared(4, 4); + ppGpEIpEI = std::make_shared>(4, 4); } void MbD::ConstVelConstraintIqcJc::calcPostDynCorrectorIteration() @@ -99,7 +99,7 @@ void MbD::ConstVelConstraintIqcJc::initialize() { ConstVelConstraintIJ::initialize(); pGpEI = std::make_shared>(4); - ppGpEIpEI = std::make_shared(4, 4); + ppGpEIpEI = std::make_shared>(4, 4); } void MbD::ConstVelConstraintIqcJc::useEquationNumbers() diff --git a/OndselSolver/ConstVelConstraintIqcJqc.cpp b/OndselSolver/ConstVelConstraintIqcJqc.cpp index 052f813..0ea8950 100644 --- a/OndselSolver/ConstVelConstraintIqcJqc.cpp +++ b/OndselSolver/ConstVelConstraintIqcJqc.cpp @@ -16,8 +16,8 @@ using namespace MbD; MbD::ConstVelConstraintIqcJqc::ConstVelConstraintIqcJqc(EndFrmsptr frmi, EndFrmsptr frmj) : ConstVelConstraintIqcJc(frmi, frmj) { pGpEJ = std::make_shared>(4); - ppGpEIpEJ = std::make_shared(4, 4); - ppGpEJpEJ = std::make_shared(4, 4); + ppGpEIpEJ = std::make_shared>(4, 4); + ppGpEJpEJ = std::make_shared>(4, 4); } void MbD::ConstVelConstraintIqcJqc::calcPostDynCorrectorIteration() @@ -121,8 +121,8 @@ void MbD::ConstVelConstraintIqcJqc::initialize() { ConstVelConstraintIqcJc::initialize(); pGpEJ = std::make_shared>(4); - ppGpEIpEJ = std::make_shared(4, 4); - ppGpEJpEJ = std::make_shared(4, 4); + ppGpEIpEJ = std::make_shared>(4, 4); + ppGpEJpEJ = std::make_shared>(4, 4); } void MbD::ConstVelConstraintIqcJqc::useEquationNumbers() diff --git a/OndselSolver/DiagonalMatrix.cpp b/OndselSolver/DiagonalMatrix.cpp index c03cfc9..13bba29 100644 --- a/OndselSolver/DiagonalMatrix.cpp +++ b/OndselSolver/DiagonalMatrix.cpp @@ -7,90 +7,5 @@ ***************************************************************************/ #include "DiagonalMatrix.h" -#include "FullMatrix.h" -namespace MbD { - - DiagMatDsptr DiagonalMatrix::times(double factor) - { - auto nrow = (int)this->size(); - auto answer = std::make_shared(nrow); - for (int i = 0; i < nrow; i++) - { - answer->at(i) = this->at(i) * factor; - } - return answer; - } - void DiagonalMatrix::atiputDiagonalMatrix(int i, std::shared_ptr diagMat) - { - for (int ii = 0; ii < diagMat->size(); ii++) - { - this->at(i + ii) = diagMat->at(ii); - } - } - FColsptr DiagonalMatrix::timesFullColumn(FColsptr fullCol) - { - //"a*b = a(i,j)b(j) sum j." - - auto nrow = (int)this->size(); - auto answer = std::make_shared>(nrow); - for (int i = 0; i < nrow; i++) - { - answer->at(i) = this->at(i) * fullCol->at(i); - } - return answer; - } - FMatDsptr DiagonalMatrix::timesFullMatrix(FMatDsptr fullMat) - { - auto nrow = (int)this->size(); - auto answer = std::make_shared(nrow); - for (int i = 0; i < nrow; i++) - { - answer->at(i) = fullMat->at(i)->times(this->at(i)); - } - return answer; - } - double DiagonalMatrix::sumOfSquares() - { - double sum = 0.0; - for (int i = 0; i < this->size(); i++) - { - double element = this->at(i); - sum += element * element; - } - return sum; - } - int DiagonalMatrix::numberOfElements() - { - auto n = (int)this->size(); - return n * n; - } - void DiagonalMatrix::zeroSelf() - { - for (int i = 0; i < this->size(); i++) { - this->at(i) = 0.0; - } - } - double DiagonalMatrix::maxMagnitude() - { - double max = 0.0; - for (int i = 0; i < this->size(); i++) - { - double element = this->at(i); - if (element < 0.0) element = -element; - if (max < element) max = element; - } - return max; - } - std::ostream& DiagonalMatrix::printOn(std::ostream& s) const - { - s << "DiagMat["; - s << this->at(0); - for (int i = 1; i < this->size(); i++) - { - s << ", " << this->at(i); - } - s << "]"; - return s; - } -} +using namespace MbD; diff --git a/OndselSolver/DiagonalMatrix.h b/OndselSolver/DiagonalMatrix.h index e802b0c..69c0c28 100644 --- a/OndselSolver/DiagonalMatrix.h +++ b/OndselSolver/DiagonalMatrix.h @@ -8,26 +8,30 @@ #pragma once -#include "FullColumn.ref.h" -#include "FullRow.ref.h" -#include "FullMatrix.ref.h" -#include "DiagonalMatrix.ref.h" #include "Array.h" #include "FullColumn.h" +#include "FullMatrix.h" namespace MbD { - class DiagonalMatrix : public Array + template + class DiagonalMatrix; + template + using DiagMatsptr = std::shared_ptr>; + using DiagMatDsptr = std::shared_ptr>; + + template + class DiagonalMatrix : public Array { // public: - DiagonalMatrix() : Array() {} - explicit DiagonalMatrix(int count) : Array(count) {} - DiagonalMatrix(int count, const double& value) : Array(count, value) {} - DiagonalMatrix(std::initializer_list list) : Array{ list } {} - void atiputDiagonalMatrix(int i, std::shared_ptr diagMat); - DiagMatDsptr times(double factor); - FColsptr timesFullColumn(FColsptr fullCol); - FMatDsptr timesFullMatrix(FMatDsptr fullMat); + DiagonalMatrix() : Array() {} + DiagonalMatrix(int count) : Array(count) {} + DiagonalMatrix(int count, const T& value) : Array(count, value) {} + DiagonalMatrix(std::initializer_list list) : Array{ list } {} + void atiputDiagonalMatrix(int i, std::shared_ptr> diagMat); + DiagMatsptr times(T factor); + FColsptr timesFullColumn(FColsptr fullCol); + FMatsptr timesFullMatrix(FMatsptr fullMat); int nrow() { return (int)this->size(); } @@ -40,6 +44,109 @@ namespace MbD { double maxMagnitude() override; std::ostream& printOn(std::ostream& s) const override; + }; + template<> + inline DiagMatDsptr DiagonalMatrix::times(double factor) + { + auto nrow = (int)this->size(); + auto answer = std::make_shared>(nrow); + for (int i = 0; i < nrow; i++) + { + answer->at(i) = this->at(i) * factor; + } + return answer; + } + template + inline void DiagonalMatrix::atiputDiagonalMatrix(int i, std::shared_ptr> diagMat) + { + for (int ii = 0; ii < diagMat->size(); ii++) + { + this->at(i + ii) = diagMat->at(ii); + } + } + template + inline DiagMatsptr DiagonalMatrix::times(T factor) + { + assert(false); + } + template + inline FColsptr DiagonalMatrix::timesFullColumn(FColsptr fullCol) + { + //"a*b = a(i,j)b(j) sum j." + + auto nrow = (int)this->size(); + auto answer = std::make_shared>(nrow); + for (int i = 0; i < nrow; i++) + { + answer->at(i) = this->at(i) * fullCol->at(i); + } + return answer; + } + template + inline FMatsptr DiagonalMatrix::timesFullMatrix(FMatsptr fullMat) + { + auto nrow = (int)this->size(); + auto answer = std::make_shared>(nrow); + for (int i = 0; i < nrow; i++) + { + answer->at(i) = fullMat->at(i)->times(this->at(i)); + } + return answer; + } + template<> + inline double DiagonalMatrix::sumOfSquares() + { + double sum = 0.0; + for (int i = 0; i < this->size(); i++) + { + double element = this->at(i); + sum += element * element; + } + return sum; + } + template + inline int DiagonalMatrix::numberOfElements() + { + auto n = (int)this->size(); + return n * n; + } + template<> + inline void DiagonalMatrix::zeroSelf() + { + for (int i = 0; i < this->size(); i++) { + this->at(i) = 0.0; + } + } + template<> + inline double DiagonalMatrix::maxMagnitude() + { + double max = 0.0; + for (int i = 0; i < this->size(); i++) + { + double element = this->at(i); + if (element < 0.0) element = -element; + if (max < element) max = element; + } + return max; + } + template + inline double DiagonalMatrix::maxMagnitude() + { + assert(false); + return 0.0; + } + template + inline std::ostream& DiagonalMatrix::printOn(std::ostream& s) const + { + s << "DiagMat["; + s << this->at(0); + for (int i = 1; i < this->size(); i++) + { + s << ", " << this->at(i); + } + s << "]"; + return s; + } } diff --git a/OndselSolver/DiagonalMatrix.ref.h b/OndselSolver/DiagonalMatrix.ref.h deleted file mode 100644 index 15127de..0000000 --- a/OndselSolver/DiagonalMatrix.ref.h +++ /dev/null @@ -1,9 +0,0 @@ -#pragma once - -#include - -namespace MbD { - class DiagonalMatrix; - using DiagMatsptr = std::shared_ptr; - using DiagMatDsptr = std::shared_ptr; -} \ No newline at end of file diff --git a/OndselSolver/DifferenceOperator.cpp b/OndselSolver/DifferenceOperator.cpp index 510a497..d71efe9 100644 --- a/OndselSolver/DifferenceOperator.cpp +++ b/OndselSolver/DifferenceOperator.cpp @@ -63,7 +63,7 @@ void DifferenceOperator::setorder(int o) void DifferenceOperator::instantiateTaylorMatrix() { if (taylorMatrix == nullptr || (taylorMatrix->nrow() != (order + 1))) { - taylorMatrix = std::make_shared(order + 1, order + 1); + taylorMatrix = std::make_shared>(order + 1, order + 1); } } diff --git a/OndselSolver/DirectionCosineIecJec.cpp b/OndselSolver/DirectionCosineIecJec.cpp index b82f874..942570f 100644 --- a/OndselSolver/DirectionCosineIecJec.cpp +++ b/OndselSolver/DirectionCosineIecJec.cpp @@ -25,7 +25,7 @@ namespace MbD { { aAjOIe = frmI->aAjOe(axisI); aAjOJe = frmJ->aAjOe(axisJ); - aAijIeJe = aAjOIe->dotVec(aAjOJe); + aAijIeJe = aAjOIe->dot(aAjOJe); } double MbD::DirectionCosineIecJec::value() diff --git a/OndselSolver/DirectionCosineIeqcJec.cpp b/OndselSolver/DirectionCosineIeqcJec.cpp index 80f438a..6e4423b 100644 --- a/OndselSolver/DirectionCosineIeqcJec.cpp +++ b/OndselSolver/DirectionCosineIeqcJec.cpp @@ -24,7 +24,7 @@ void DirectionCosineIeqcJec::initialize() { DirectionCosineIecJec::initialize(); pAijIeJepEI = std::make_shared>(4); - ppAijIeJepEIpEI = std::make_shared(4, 4); + ppAijIeJepEIpEI = std::make_shared>(4, 4); } void DirectionCosineIeqcJec::initializeGlobally() @@ -56,7 +56,7 @@ void DirectionCosineIeqcJec::calcPostDynCorrectorIteration() auto& ppAjOIepEIipEI = ppAjOIepEIpEI->at(i); for (int j = 0; j < 4; j++) { - ppAijIeJepEIipEI->at(j) = ppAjOIepEIipEI->at(j)->dotVec(aAjOJe); + ppAijIeJepEIipEI->at(j) = ppAjOIepEIipEI->at(j)->dot(aAjOJe); } } ppAijIeJepEIpEI->symLowerWithUpper(); diff --git a/OndselSolver/DirectionCosineIeqcJeqc.cpp b/OndselSolver/DirectionCosineIeqcJeqc.cpp index e449642..92598c9 100644 --- a/OndselSolver/DirectionCosineIeqcJeqc.cpp +++ b/OndselSolver/DirectionCosineIeqcJeqc.cpp @@ -24,8 +24,8 @@ void DirectionCosineIeqcJeqc::initialize() { DirectionCosineIeqcJec::initialize(); pAijIeJepEJ = std::make_shared>(4); - ppAijIeJepEIpEJ = std::make_shared(4, 4); - ppAijIeJepEJpEJ = std::make_shared(4, 4); + ppAijIeJepEIpEJ = std::make_shared>(4, 4); + ppAijIeJepEJpEJ = std::make_shared>(4, 4); } void DirectionCosineIeqcJeqc::initializeGlobally() @@ -45,7 +45,7 @@ void DirectionCosineIeqcJeqc::calcPostDynCorrectorIteration() pAjOJepEJT = std::static_pointer_cast(frmJ)->pAjOepET(axisJ); for (int i = 0; i < 4; i++) { - pAijIeJepEJ->at(i) = aAjOIe->dotVec(pAjOJepEJT->at(i)); + pAijIeJepEJ->at(i) = aAjOIe->dot(pAjOJepEJT->at(i)); } for (int i = 0; i < 4; i++) { @@ -61,7 +61,7 @@ void DirectionCosineIeqcJeqc::calcPostDynCorrectorIteration() auto& ppAjOJepEJipEJ = ppAjOJepEJpEJ->at(i); for (int j = 0; j < 4; j++) { - ppAijIeJepEJipEJ->at(j) = aAjOIe->dotVec(ppAjOJepEJipEJ->at(j)); + ppAijIeJepEJipEJ->at(j) = aAjOIe->dot(ppAjOJepEJipEJ->at(j)); } } ppAijIeJepEJpEJ->symLowerWithUpper(); diff --git a/OndselSolver/DirectionCosineIeqctJeqc.cpp b/OndselSolver/DirectionCosineIeqctJeqc.cpp index 5790b55..6fb505c 100644 --- a/OndselSolver/DirectionCosineIeqctJeqc.cpp +++ b/OndselSolver/DirectionCosineIeqctJeqc.cpp @@ -60,7 +60,7 @@ void DirectionCosineIeqctJeqc::preVelIC() { Item::preVelIC(); auto pAjOIept = std::static_pointer_cast(frmI)->pAjOept(axisI); - pAijIeJept = pAjOIept->dotVec(aAjOJe); + pAijIeJept = pAjOIept->dot(aAjOJe); } double DirectionCosineIeqctJeqc::pvaluept() @@ -82,7 +82,7 @@ void DirectionCosineIeqctJeqc::preAccIC() } for (int i = 0; i < 4; i++) { - ppAijIeJepEJpt->atiput(i, pAjOIept->dotVec(pAjOJepEJT->at(i))); + ppAijIeJepEJpt->atiput(i, pAjOIept->dot(pAjOJepEJT->at(i))); } - ppAijIeJeptpt = ppAjOIeptpt->dotVec(aAjOJe); + ppAijIeJeptpt = ppAjOIeptpt->dot(aAjOJe); } diff --git a/OndselSolver/DispCompIecJecIe.cpp b/OndselSolver/DispCompIecJecIe.cpp index 515a6fa..213a0ff 100644 --- a/OndselSolver/DispCompIecJecIe.cpp +++ b/OndselSolver/DispCompIecJecIe.cpp @@ -23,7 +23,7 @@ void MbD::DispCompIecJecIe::calc_value() { aAjOIe = frmI->aAjOe(axis); rIeJeO = frmJ->rOeO->minusFullColumn(frmI->rOeO); - riIeJeIe = aAjOIe->dotVec(rIeJeO); + riIeJeIe = aAjOIe->dot(rIeJeO); } void MbD::DispCompIecJecIe::calcPostDynCorrectorIteration() diff --git a/OndselSolver/DispCompIecJecKeqc.cpp b/OndselSolver/DispCompIecJecKeqc.cpp index 7584dc6..efe61f5 100644 --- a/OndselSolver/DispCompIecJecKeqc.cpp +++ b/OndselSolver/DispCompIecJecKeqc.cpp @@ -22,7 +22,7 @@ DispCompIecJecKeqc::DispCompIecJecKeqc(EndFrmsptr frmi, EndFrmsptr frmj, EndFrms void DispCompIecJecKeqc::initialize() { priIeJeKepEK = std::make_shared>(4); - ppriIeJeKepEKpEK = std::make_shared(4, 4); + ppriIeJeKepEKpEK = std::make_shared>(4, 4); } void DispCompIecJecKeqc::initializeGlobally() @@ -37,7 +37,7 @@ void DispCompIecJecKeqc::calcPostDynCorrectorIteration() auto efrmKqc = std::static_pointer_cast(efrmK); aAjOKe = efrmKqc->aAjOe(axisK); rIeJeO = frmJqc->rOeO->minusFullColumn(frmIqc->rOeO); - riIeJeKe = aAjOKe->dotVec(rIeJeO); + riIeJeKe = aAjOKe->dot(rIeJeO); pAjOKepEKT = efrmKqc->pAjOepET(axisK); ppAjOKepEKpEK = efrmKqc->ppAjOepEpE(axisK); for (int i = 0; i < 4; i++) @@ -45,10 +45,10 @@ void DispCompIecJecKeqc::calcPostDynCorrectorIteration() priIeJeKepEK->at(i) = ((pAjOKepEKT->at(i))->dot(rIeJeO)); auto& ppAjOKepEKipEK = ppAjOKepEKpEK->at(i); auto& ppriIeJeKepEKipEK = ppriIeJeKepEKpEK->at(i); - ppriIeJeKepEKipEK->at(i) = ((ppAjOKepEKipEK->at(i))->dotVec(rIeJeO)); + ppriIeJeKepEKipEK->at(i) = ((ppAjOKepEKipEK->at(i))->dot(rIeJeO)); for (int j = i + 1; j < 4; j++) { - auto ppriIeJeKepEKipEKj = (ppAjOKepEKipEK->at(i))->dotVec(rIeJeO); + auto ppriIeJeKepEKipEKj = (ppAjOKepEKipEK->at(i))->dot(rIeJeO); ppriIeJeKepEKipEK->at(j) = ppriIeJeKepEKipEKj; ppriIeJeKepEKpEK->at(j)->at(i) = ppriIeJeKepEKipEKj; } diff --git a/OndselSolver/DispCompIeqcJecIe.cpp b/OndselSolver/DispCompIeqcJecIe.cpp index ec8c065..8a684a7 100644 --- a/OndselSolver/DispCompIeqcJecIe.cpp +++ b/OndselSolver/DispCompIeqcJecIe.cpp @@ -19,8 +19,8 @@ MbD::DispCompIeqcJecIe::DispCompIeqcJecIe(EndFrmsptr frmi, EndFrmsptr frmj, int { priIeJeIepXI = std::make_shared>(3); priIeJeIepEI = std::make_shared>(4); - ppriIeJeIepXIpEI = std::make_shared(3, 4); - ppriIeJeIepEIpEI = std::make_shared(4, 4); + ppriIeJeIepXIpEI = std::make_shared>(3, 4); + ppriIeJeIepEIpEI = std::make_shared>(4, 4); } void MbD::DispCompIeqcJecIe::calc_ppvaluepEIpEI() @@ -35,10 +35,10 @@ void MbD::DispCompIeqcJecIe::calc_ppvaluepEIpEI() auto ppriIeJeIepEIipEI = ppriIeJeIepEIpEI->at(i); for (int j = i; j < 4; j++) { - auto term1 = ppAjOIepEIipEI->at(j)->dotVec(rIeJeO); + auto term1 = ppAjOIepEIipEI->at(j)->dot(rIeJeO); auto mterm2 = pAjOIepEIT->at(i)->dot(mprIeJeOpEIT->at(j)); auto mterm3 = (i == j) ? mterm2 : pAjOIepEIT->at(j)->dot(mprIeJeOpEIT->at(i)); - auto mterm4 = aAjOIe->dotVec(mpprIeJeOpEIipEI->at(j)); + auto mterm4 = aAjOIe->dot(mpprIeJeOpEIipEI->at(j)); ppriIeJeIepEIipEI->atiput(j, term1 - mterm2 - mterm3 - mterm4); } } @@ -64,7 +64,7 @@ void MbD::DispCompIeqcJecIe::calc_pvaluepEI() auto mprIeJeOpEIT = frmIeqc->prOeOpE->transpose(); for (int i = 0; i < 4; i++) { - priIeJeIepEI->atiput(i, pAjOIepEIT->at(i)->dot(rIeJeO) - aAjOIe->dotVec(mprIeJeOpEIT->at(i))); + priIeJeIepEI->atiput(i, pAjOIepEIT->at(i)->dot(rIeJeO) - aAjOIe->dot(mprIeJeOpEIT->at(i))); } } @@ -91,8 +91,8 @@ void MbD::DispCompIeqcJecIe::initialize() DispCompIecJecIe::initialize(); priIeJeIepXI = std::make_shared>(3); priIeJeIepEI = std::make_shared>(4); - ppriIeJeIepXIpEI = std::make_shared(3, 4); - ppriIeJeIepEIpEI = std::make_shared(4, 4); + ppriIeJeIepXIpEI = std::make_shared>(3, 4); + ppriIeJeIepEIpEI = std::make_shared>(4, 4); } void MbD::DispCompIeqcJecIe::initializeGlobally() diff --git a/OndselSolver/DispCompIeqcJecKeqc.cpp b/OndselSolver/DispCompIeqcJecKeqc.cpp index b83fd54..c25efc1 100644 --- a/OndselSolver/DispCompIeqcJecKeqc.cpp +++ b/OndselSolver/DispCompIeqcJecKeqc.cpp @@ -24,9 +24,9 @@ void DispCompIeqcJecKeqc::initialize() DispCompIecJecKeqc::initialize(); priIeJeKepXI = std::make_shared>(3); priIeJeKepEI = std::make_shared>(4); - ppriIeJeKepEIpEI = std::make_shared(4, 4); - ppriIeJeKepXIpEK = std::make_shared(3, 4); - ppriIeJeKepEIpEK = std::make_shared(4, 4); + ppriIeJeKepEIpEI = std::make_shared>(4, 4); + ppriIeJeKepXIpEK = std::make_shared>(3, 4); + ppriIeJeKepEIpEK = std::make_shared>(4, 4); } void DispCompIeqcJecKeqc::calcPostDynCorrectorIteration() @@ -41,7 +41,7 @@ void DispCompIeqcJecKeqc::calcPostDynCorrectorIteration() } for (int i = 0; i < 4; i++) { - priIeJeKepEI->at(i) = 0.0 - (aAjOKe->dotVec(mprIeJeOpEIT->at(i))); + priIeJeKepEI->at(i) = 0.0 - (aAjOKe->dot(mprIeJeOpEIT->at(i))); } for (int i = 0; i < 3; i++) { @@ -55,10 +55,10 @@ void DispCompIeqcJecKeqc::calcPostDynCorrectorIteration() { auto& mpprIeJeOpEIipEI = mpprIeJeOpEIpEI->at(i); auto& ppriIeJeKepEIipEI = ppriIeJeKepEIpEI->at(i); - ppriIeJeKepEIipEI->at(i) = 0.0 - (aAjOKe->dotVec(mpprIeJeOpEIipEI->at(i))); + ppriIeJeKepEIipEI->at(i) = 0.0 - (aAjOKe->dot(mpprIeJeOpEIipEI->at(i))); for (int j = 0; j < 4; j++) { - auto ppriIeJeKepEIipEIj = 0.0 - (aAjOKe->dotVec(mpprIeJeOpEIipEI->at(j))); + auto ppriIeJeKepEIipEIj = 0.0 - (aAjOKe->dot(mpprIeJeOpEIipEI->at(j))); ppriIeJeKepEIipEI->at(j) = ppriIeJeKepEIipEIj; ppriIeJeKepEIpEI->at(j)->at(i) = ppriIeJeKepEIipEIj; } diff --git a/OndselSolver/DispCompIeqcJeqcIe.cpp b/OndselSolver/DispCompIeqcJeqcIe.cpp index 33a40bf..662b43e 100644 --- a/OndselSolver/DispCompIeqcJeqcIe.cpp +++ b/OndselSolver/DispCompIeqcJeqcIe.cpp @@ -19,9 +19,9 @@ MbD::DispCompIeqcJeqcIe::DispCompIeqcJeqcIe(EndFrmsptr frmi, EndFrmsptr frmj, in { priIeJeIepXJ = std::make_shared>(3); priIeJeIepEJ = std::make_shared>(4); - ppriIeJeIepEIpXJ = std::make_shared(4, 3); - ppriIeJeIepEIpEJ = std::make_shared(4, 4); - ppriIeJeIepEJpEJ = std::make_shared(4, 4); + ppriIeJeIepEIpXJ = std::make_shared>(4, 3); + ppriIeJeIepEIpEJ = std::make_shared>(4, 4); + ppriIeJeIepEJpEJ = std::make_shared>(4, 4); } void MbD::DispCompIeqcJeqcIe::calc_ppvaluepEIpEJ() @@ -46,7 +46,7 @@ void MbD::DispCompIeqcJeqcIe::calc_ppvaluepEJpEJ() auto ppriIeJeIepEJipEJ = ppriIeJeIepEJpEJ->at(i); for (int j = i; j < 4; j++) { - auto term1 = aAjOIe->dotVec(pprIeJeOpEJipEJ->at(j)); + auto term1 = aAjOIe->dot(pprIeJeOpEJipEJ->at(j)); ppriIeJeIepEJipEJ->atiput(j, term1); } } @@ -59,7 +59,7 @@ void MbD::DispCompIeqcJeqcIe::calc_pvaluepEJ() auto prIeJeOpEJT = frmJeqc->prOeOpE->transpose(); for (int i = 0; i < 4; i++) { - priIeJeIepEJ->atiput(i, aAjOIe->dotVec(prIeJeOpEJT->at(i))); + priIeJeIepEJ->atiput(i, aAjOIe->dot(prIeJeOpEJT->at(i))); } } @@ -88,9 +88,9 @@ void MbD::DispCompIeqcJeqcIe::initialize() DispCompIeqcJecIe::initialize(); priIeJeIepXJ = std::make_shared>(3); priIeJeIepEJ = std::make_shared>(4); - ppriIeJeIepEIpXJ = std::make_shared(4, 3); - ppriIeJeIepEIpEJ = std::make_shared(4, 4); - ppriIeJeIepEJpEJ = std::make_shared(4, 4); + ppriIeJeIepEIpXJ = std::make_shared>(4, 3); + ppriIeJeIepEIpEJ = std::make_shared>(4, 4); + ppriIeJeIepEJpEJ = std::make_shared>(4, 4); } FMatDsptr MbD::DispCompIeqcJeqcIe::ppvaluepEIpEJ() diff --git a/OndselSolver/DispCompIeqcJeqcKeqc.cpp b/OndselSolver/DispCompIeqcJeqcKeqc.cpp index 72210be..61a6254 100644 --- a/OndselSolver/DispCompIeqcJeqcKeqc.cpp +++ b/OndselSolver/DispCompIeqcJeqcKeqc.cpp @@ -24,9 +24,9 @@ void DispCompIeqcJeqcKeqc::initialize() DispCompIeqcJecKeqc::initialize(); priIeJeKepXJ = std::make_shared>(3); priIeJeKepEJ = std::make_shared>(4); - ppriIeJeKepEJpEJ = std::make_shared(4, 4); - ppriIeJeKepXJpEK = std::make_shared(3, 4); - ppriIeJeKepEJpEK = std::make_shared(4, 4); + ppriIeJeKepEJpEJ = std::make_shared>(4, 4); + ppriIeJeKepXJpEK = std::make_shared>(3, 4); + ppriIeJeKepEJpEK = std::make_shared>(4, 4); } void DispCompIeqcJeqcKeqc::calcPostDynCorrectorIteration() @@ -41,7 +41,7 @@ void DispCompIeqcJeqcKeqc::calcPostDynCorrectorIteration() } for (int i = 0; i < 4; i++) { - priIeJeKepEJ->atiput(i, aAjOKe->dotVec(prIeJeOpEJT->at(i))); + priIeJeKepEJ->atiput(i, aAjOKe->dot(prIeJeOpEJT->at(i))); } for (int i = 0; i < 3; i++) { @@ -55,10 +55,10 @@ void DispCompIeqcJeqcKeqc::calcPostDynCorrectorIteration() { auto& pprIeJeOpEJipEJ = pprIeJeOpEJpEJ->at(i); auto& ppriIeJeKepEJipEJ = ppriIeJeKepEJpEJ->at(i); - ppriIeJeKepEJipEJ->atiput(i, aAjOKe->dotVec(pprIeJeOpEJipEJ->at(i))); + ppriIeJeKepEJipEJ->atiput(i, aAjOKe->dot(pprIeJeOpEJipEJ->at(i))); for (int j = 0; j < 4; j++) { - auto ppriIeJeKepEJipEJj = (aAjOKe->dotVec(pprIeJeOpEJipEJ->at(j))); + auto ppriIeJeKepEJipEJj = (aAjOKe->dot(pprIeJeOpEJipEJ->at(j))); ppriIeJeKepEJipEJ->atiput(j, ppriIeJeKepEJipEJj); ppriIeJeKepEJpEJ->atijput(j, i, ppriIeJeKepEJipEJj); } diff --git a/OndselSolver/DispCompIeqcJeqcKeqct.cpp b/OndselSolver/DispCompIeqcJeqcKeqct.cpp index 3a1b31a..0820af1 100644 --- a/OndselSolver/DispCompIeqcJeqcKeqct.cpp +++ b/OndselSolver/DispCompIeqcJeqcKeqct.cpp @@ -47,7 +47,7 @@ void DispCompIeqcJeqcKeqct::preVelIC() { Item::preVelIC(); auto pAjOKept = std::static_pointer_cast(efrmK)->pAjOept(axisK); - priIeJeKept = pAjOKept->dotVec(rIeJeO); + priIeJeKept = pAjOKept->dot(rIeJeO); } double DispCompIeqcJeqcKeqct::pvaluept() @@ -100,9 +100,9 @@ void DispCompIeqcJeqcKeqct::preAccIC() } for (int i = 0; i < 4; i++) { - ppriIeJeKepEIpt->atiput(i, pAjOKept->dotVec(prIeJeOpEIT->at(i))); - ppriIeJeKepEJpt->atiput(i, pAjOKept->dotVec(prIeJeOpEJT->at(i))); + ppriIeJeKepEIpt->atiput(i, pAjOKept->dot(prIeJeOpEIT->at(i))); + ppriIeJeKepEJpt->atiput(i, pAjOKept->dot(prIeJeOpEJT->at(i))); ppriIeJeKepEKpt->atiput(i, ppAjOKepEKTpt->at(i)->dot(rIeJeO)); } - ppriIeJeKeptpt = ppAjOKeptpt->dotVec(rIeJeO); + ppriIeJeKeptpt = ppAjOKeptpt->dot(rIeJeO); } diff --git a/OndselSolver/DispCompIeqctJeqcIe.cpp b/OndselSolver/DispCompIeqctJeqcIe.cpp index a4b8a9c..4b6b20d 100644 --- a/OndselSolver/DispCompIeqctJeqcIe.cpp +++ b/OndselSolver/DispCompIeqctJeqcIe.cpp @@ -38,8 +38,8 @@ void MbD::DispCompIeqctJeqcIe::calc_ppvaluepEIpt() auto mprIeJeOpEITi = mprIeJeOpEIT->at(i); auto mpprIeJeOpEITpti = mpprIeJeOpEITpt->at(i); auto ppriIeJeIepEIpti = ppAjOIepEITpti->dot(rIeJeO) - pAjOIepEITi->dot(mprIeJeOpt) - - pAjOIept->dotVec(mprIeJeOpEITi) - - aAjOIe->dotVec(mpprIeJeOpEITpti); + pAjOIept->dot(mprIeJeOpEITi) - + aAjOIe->dot(mpprIeJeOpEITpti); ppriIeJeIepEIpt->atiput(i, ppriIeJeIepEIpti); } } @@ -52,7 +52,7 @@ void MbD::DispCompIeqctJeqcIe::calc_ppvaluepEJpt() auto prIeJeOpEJT = frmJeqct->prOeOpE->transpose(); for (int i = 0; i < 4; i++) { - ppriIeJeIepEJpt->atiput(i, pAjOIept->dotVec(prIeJeOpEJT->at(i))); + ppriIeJeIepEJpt->atiput(i, pAjOIept->dot(prIeJeOpEJT->at(i))); } } @@ -63,8 +63,8 @@ void MbD::DispCompIeqctJeqcIe::calc_ppvalueptpt() auto ppAjOIeptpt = frmIeqct->ppAjOeptpt(axis); auto mprIeJeOpt = frmIeqct->prOeOpt; auto mpprIeJeOptpt = frmIeqct->pprOeOptpt; - ppriIeJeIeptpt = ppAjOIeptpt->dotVec(rIeJeO) - pAjOIept->dotVec(mprIeJeOpt) - pAjOIept->dotVec(mprIeJeOpt) - - aAjOIe->dotVec(mpprIeJeOptpt); + ppriIeJeIeptpt = ppAjOIeptpt->dot(rIeJeO) - pAjOIept->dot(mprIeJeOpt) - pAjOIept->dot(mprIeJeOpt) - + aAjOIe->dot(mpprIeJeOptpt); } void MbD::DispCompIeqctJeqcIe::calc_ppvaluepXIpt() @@ -92,7 +92,7 @@ void MbD::DispCompIeqctJeqcIe::calc_pvaluept() auto frmIeqct = std::static_pointer_cast(frmI); auto pAjOIept = frmIeqct->pAjOept(axis); auto mprIeJeOpt = frmIeqct->prOeOpt; - priIeJeIept = pAjOIept->dotVec(rIeJeO) - aAjOIe->dotVec(mprIeJeOpt); + priIeJeIept = pAjOIept->dot(rIeJeO) - aAjOIe->dot(mprIeJeOpt); } void MbD::DispCompIeqctJeqcIe::calcPostDynCorrectorIteration() diff --git a/OndselSolver/DispCompIeqctJeqcKeqct.cpp b/OndselSolver/DispCompIeqctJeqcKeqct.cpp index e3013e6..41ea675 100644 --- a/OndselSolver/DispCompIeqctJeqcKeqct.cpp +++ b/OndselSolver/DispCompIeqctJeqcKeqct.cpp @@ -23,7 +23,7 @@ void DispCompIeqctJeqcKeqct::preVelIC() { DispCompIeqcJeqcKeqct::preVelIC(); auto& mprIeJeOpt = std::static_pointer_cast(frmI)->prOeOpt; - priIeJeKept -= aAjOKe->dotVec(mprIeJeOpt); + priIeJeKept -= aAjOKe->dot(mprIeJeOpt); } void DispCompIeqctJeqcKeqct::preAccIC() @@ -36,8 +36,8 @@ void DispCompIeqctJeqcKeqct::preAccIC() auto& mpprIeJeOptpt = efrmIqct->pprOeOptpt; for (int i = 0; i < 4; i++) { - ppriIeJeKepEIpt->atiminusNumber(i, aAjOKe->dotVec(mpprIeJeOpEITpt->at(i))); + ppriIeJeKepEIpt->atiminusNumber(i, aAjOKe->dot(mpprIeJeOpEITpt->at(i))); ppriIeJeKepEKpt->atiminusNumber(i, pAjOKepEKT->at(i)->dot(mprIeJeOpt)); } - ppriIeJeKeptpt += -(2.0 * pAjOKept->dotVec(mprIeJeOpt)) - aAjOKe->dotVec(mpprIeJeOptpt); + ppriIeJeKeptpt += -(2.0 * pAjOKept->dot(mprIeJeOpt)) - aAjOKe->dot(mpprIeJeOptpt); } diff --git a/OndselSolver/DistIeqcJec.cpp b/OndselSolver/DistIeqcJec.cpp index 282b52f..9c23e35 100644 --- a/OndselSolver/DistIeqcJec.cpp +++ b/OndselSolver/DistIeqcJec.cpp @@ -41,7 +41,7 @@ void MbD::DistIeqcJec::calcPrivate() pprIeJepXIipXI->atiput(j, element / rIeJe); } } - pprIeJepXIpEI = std::make_shared(3, 4); + pprIeJepXIpEI = std::make_shared>(3, 4); for (int i = 0; i < 3; i++) { auto& pprIeJepXIipEI = pprIeJepXIpEI->at(i); @@ -53,7 +53,7 @@ void MbD::DistIeqcJec::calcPrivate() pprIeJepXIipEI->atiput(j, element / rIeJe); } } - pprIeJepEIpEI = std::make_shared(4, 4); + pprIeJepEIpEI = std::make_shared>(4, 4); for (int i = 0; i < 4; i++) { auto& pprIeJepEIipEI = pprIeJepEIpEI->at(i); @@ -63,7 +63,7 @@ void MbD::DistIeqcJec::calcPrivate() for (int j = 0; j < 4; j++) { auto element = mprIeJeOpEIiT->dot(mprIeJeOpEIT->at(j)) - - mpprIeJeOpEIipEI->at(j)->dotVec(rIeJeO) - prIeJepEIi * prIeJepEI->at(j); + - mpprIeJeOpEIipEI->at(j)->dot(rIeJeO) - prIeJepEIi * prIeJepEI->at(j); pprIeJepEIipEI->atiput(j, element / rIeJe); } } @@ -74,9 +74,9 @@ void MbD::DistIeqcJec::initialize() DistIecJec::initialize(); prIeJepXI = std::make_shared>(3); prIeJepEI = std::make_shared>(4); - pprIeJepXIpXI = std::make_shared(3, 3); - pprIeJepXIpEI = std::make_shared(3, 4); - pprIeJepEIpEI = std::make_shared(4, 4); + pprIeJepXIpXI = std::make_shared>(3, 3); + pprIeJepXIpEI = std::make_shared>(3, 4); + pprIeJepEIpEI = std::make_shared>(4, 4); } FMatDsptr MbD::DistIeqcJec::ppvaluepEIpEI() diff --git a/OndselSolver/DistIeqcJeqc.cpp b/OndselSolver/DistIeqcJeqc.cpp index d7fdb1c..8b8a3da 100644 --- a/OndselSolver/DistIeqcJeqc.cpp +++ b/OndselSolver/DistIeqcJeqc.cpp @@ -111,7 +111,7 @@ void MbD::DistIeqcJeqc::calcPrivate() for (int j = 0; j < 4; j++) { auto element = prIeJeOpEJiT->dot(prIeJeOpEJT->at(j)) - + pprIeJeOpEJipEJ->at(j)->dotVec(rIeJeO) - prIeJepEJi * prIeJepEJ->at(j); + + pprIeJeOpEJipEJ->at(j)->dot(rIeJeO) - prIeJepEJi * prIeJepEJ->at(j); pprIeJepEJipEJ->atiput(j, element / rIeJe); } } @@ -122,13 +122,13 @@ void MbD::DistIeqcJeqc::initialize() DistIeqcJec::initialize(); prIeJepXJ = std::make_shared>(3); prIeJepEJ = std::make_shared>(4); - pprIeJepXIpXJ = std::make_shared(3, 3); - pprIeJepEIpXJ = std::make_shared(4, 3); - pprIeJepXJpXJ = std::make_shared(3, 3); - pprIeJepXIpEJ = std::make_shared(3, 4); - pprIeJepEIpEJ = std::make_shared(4, 4); - pprIeJepXJpEJ = std::make_shared(3, 4); - pprIeJepEJpEJ = std::make_shared(4, 4); + pprIeJepXIpXJ = std::make_shared>(3, 3); + pprIeJepEIpXJ = std::make_shared>(4, 3); + pprIeJepXJpXJ = std::make_shared>(3, 3); + pprIeJepXIpEJ = std::make_shared>(3, 4); + pprIeJepEIpEJ = std::make_shared>(4, 4); + pprIeJepXJpEJ = std::make_shared>(3, 4); + pprIeJepEJpEJ = std::make_shared>(4, 4); } FMatDsptr MbD::DistIeqcJeqc::ppvaluepEIpEJ() diff --git a/OndselSolver/DistxyIeqcJec.cpp b/OndselSolver/DistxyIeqcJec.cpp index 165358a..12436c2 100644 --- a/OndselSolver/DistxyIeqcJec.cpp +++ b/OndselSolver/DistxyIeqcJec.cpp @@ -161,9 +161,9 @@ void MbD::DistxyIeqcJec::initialize() DistxyIecJec::initialize(); pdistxypXI = std::make_shared>(3); pdistxypEI = std::make_shared>(4); - ppdistxypXIpXI = std::make_shared(3, 3); - ppdistxypXIpEI = std::make_shared(3, 4); - ppdistxypEIpEI = std::make_shared(4, 4); + ppdistxypXIpXI = std::make_shared>(3, 3); + ppdistxypXIpEI = std::make_shared>(3, 4); + ppdistxypEIpEI = std::make_shared>(4, 4); } FMatDsptr MbD::DistxyIeqcJec::ppvaluepEIpEI() diff --git a/OndselSolver/DistxyIeqcJeqc.cpp b/OndselSolver/DistxyIeqcJeqc.cpp index f422f2d..92140ff 100644 --- a/OndselSolver/DistxyIeqcJeqc.cpp +++ b/OndselSolver/DistxyIeqcJeqc.cpp @@ -293,13 +293,13 @@ void MbD::DistxyIeqcJeqc::initialize() DistxyIeqcJec::initialize(); pdistxypXJ = std::make_shared>(3); pdistxypEJ = std::make_shared>(4); - ppdistxypXIpXJ = std::make_shared(3, 3); - ppdistxypXIpEJ = std::make_shared(3, 4); - ppdistxypEIpXJ = std::make_shared(4, 3); - ppdistxypEIpEJ = std::make_shared(4, 4); - ppdistxypXJpXJ = std::make_shared(3, 3); - ppdistxypXJpEJ = std::make_shared(3, 4); - ppdistxypEJpEJ = std::make_shared(4, 4); + ppdistxypXIpXJ = std::make_shared>(3, 3); + ppdistxypXIpEJ = std::make_shared>(3, 4); + ppdistxypEIpXJ = std::make_shared>(4, 3); + ppdistxypEIpEJ = std::make_shared>(4, 4); + ppdistxypXJpXJ = std::make_shared>(3, 3); + ppdistxypXJpEJ = std::make_shared>(3, 4); + ppdistxypEJpEJ = std::make_shared>(4, 4); } FMatDsptr MbD::DistxyIeqcJeqc::ppvaluepEIpEJ() diff --git a/OndselSolver/EndFramec.h b/OndselSolver/EndFramec.h index 9e25aba..d370048 100644 --- a/OndselSolver/EndFramec.h +++ b/OndselSolver/EndFramec.h @@ -45,7 +45,7 @@ namespace MbD { MarkerFrame* markerFrame; //Use raw pointer when pointing backwards. FColDsptr rOeO = std::make_shared>(3); - FMatDsptr aAOe = FullMatrixDouble::identitysptr(3); + FMatDsptr aAOe = FullMatrix::identitysptr(3); }; //using EndFrmsptr = std::shared_ptr; } diff --git a/OndselSolver/EndFrameqc.cpp b/OndselSolver/EndFrameqc.cpp index 746e22a..f938f12 100644 --- a/OndselSolver/EndFrameqc.cpp +++ b/OndselSolver/EndFrameqc.cpp @@ -25,10 +25,10 @@ EndFrameqc::EndFrameqc(const char* str) : EndFramec(str) { void EndFrameqc::initialize() { - prOeOpE = std::make_shared(3, 4); - pprOeOpEpE = std::make_shared(4, 4); + prOeOpE = std::make_shared>(3, 4); + pprOeOpEpE = std::make_shared>(4, 4); pAOepE = std::make_shared>(4); - ppAOepEpE = std::make_shared(4, 4); + ppAOepEpE = std::make_shared>(4, 4); } void EndFrameqc::initializeGlobally() @@ -59,7 +59,7 @@ void MbD::EndFrameqc::initEndFrameqct2() FMatFColDsptr EndFrameqc::ppAjOepEpE(int jj) { - auto answer = std::make_shared(4, 4); + auto answer = std::make_shared>(4, 4); for (int i = 0; i < 4; i++) { auto& answeri = answer->at(i); auto& ppAOepEipE = ppAOepEpE->at(i); @@ -80,7 +80,7 @@ void EndFrameqc::calcPostDynCorrectorIteration() FMatDsptr EndFrameqc::pAjOepET(int axis) { - auto answer = std::make_shared(4, 3); + auto answer = std::make_shared>(4, 3); for (int i = 0; i < 4; i++) { auto& answeri = answer->at(i); auto& pAOepEi = pAOepE->at(i); @@ -94,7 +94,7 @@ FMatDsptr EndFrameqc::pAjOepET(int axis) FMatDsptr EndFrameqc::ppriOeOpEpE(int ii) { - auto answer = std::make_shared(4, 4); + auto answer = std::make_shared>(4, 4); for (int i = 0; i < 4; i++) { auto& answeri = answer->at(i); auto& pprOeOpEipE = pprOeOpEpE->at(i); diff --git a/OndselSolver/EndFrameqct.cpp b/OndselSolver/EndFrameqct.cpp index bdd3a67..efd7b34 100644 --- a/OndselSolver/EndFrameqct.cpp +++ b/OndselSolver/EndFrameqct.cpp @@ -9,6 +9,7 @@ #include "EndFrameqct.h" #include "MarkerFrame.h" #include "System.h" +#include "Symbolic.h" #include "Time.h" #include "EulerParameters.h" #include "CREATE.h" @@ -16,312 +17,352 @@ #include "EulerAngleszxzDot.h" #include "EulerAngleszxzDDot.h" -namespace MbD { - template class EulerParameters; +using namespace MbD; - EndFrameqct::EndFrameqct() { - } - - EndFrameqct::EndFrameqct(const char *str) : EndFrameqc(str) { - } - - void EndFrameqct::initialize() { - EndFrameqc::initialize(); - rmem = std::make_shared>(3); - prmempt = std::make_shared>(3); - pprmemptpt = std::make_shared>(3); - aAme = FullMatrixDouble::identitysptr(3); - pAmept = std::make_shared(3, 3); - ppAmeptpt = std::make_shared(3, 3); - pprOeOpEpt = std::make_shared(3, 4); - pprOeOptpt = std::make_shared>(3); - ppAOepEpt = std::make_shared>(4); - ppAOeptpt = std::make_shared(3, 3); - } - - void EndFrameqct::initializeLocally() { - if (!rmemBlks) { - rmem->zeroSelf(); - prmempt->zeroSelf(); - pprmemptpt->zeroSelf(); - } - if (!phiThePsiBlks) { - aAme->identity(); - pAmept->zeroSelf(); - ppAmeptpt->zeroSelf(); - } - } - - void EndFrameqct::initializeGlobally() { - if (rmemBlks) { - initprmemptBlks(); - initpprmemptptBlks(); - } - if (phiThePsiBlks) { - initpPhiThePsiptBlks(); - initppPhiThePsiptptBlks(); - } - } - - void EndFrameqct::initprmemptBlks() { - auto &mbdTime = this->root()->time; - prmemptBlks = std::make_shared>(3); - for (int i = 0; i < 3; i++) { - auto &disp = rmemBlks->at(i); - auto var = disp->differentiateWRT(mbdTime); - auto vel = var->simplified(var); - prmemptBlks->at(i) = vel; - } - } - - void EndFrameqct::initpprmemptptBlks() { - auto &mbdTime = this->root()->time; - pprmemptptBlks = std::make_shared>(3); - for (int i = 0; i < 3; i++) { - auto &vel = prmemptBlks->at(i); - auto var = vel->differentiateWRT(mbdTime); - auto acc = var->simplified(var); - pprmemptptBlks->at(i) = acc; - } - } - - void EndFrameqct::initpPhiThePsiptBlks() { - auto &mbdTime = this->root()->time; - pPhiThePsiptBlks = std::make_shared>(3); - for (int i = 0; i < 3; i++) { - auto &angle = phiThePsiBlks->at(i); - auto var = angle->differentiateWRT(mbdTime); - //std::cout << "var " << *var << std::endl; - auto vel = var->simplified(var); - //std::cout << "vel " << *vel << std::endl; - pPhiThePsiptBlks->at(i) = vel; - } - } - - void EndFrameqct::initppPhiThePsiptptBlks() { - auto &mbdTime = this->root()->time; - ppPhiThePsiptptBlks = std::make_shared>(3); - for (int i = 0; i < 3; i++) { - auto &angleVel = pPhiThePsiptBlks->at(i); - auto var = angleVel->differentiateWRT(mbdTime); - auto angleAcc = var->simplified(var); - ppPhiThePsiptptBlks->at(i) = angleAcc; - } - } - - void EndFrameqct::postInput() { - this->evalrmem(); - this->evalAme(); - Item::postInput(); - } - - void EndFrameqct::calcPostDynCorrectorIteration() { - auto &rOmO = markerFrame->rOmO; - auto &aAOm = markerFrame->aAOm; - rOeO = rOmO->plusFullColumn(aAOm->timesFullColumn(rmem)); - auto &prOmOpE = markerFrame->prOmOpE; - auto &pAOmpE = markerFrame->pAOmpE; - for (int i = 0; i < 3; i++) { - auto &prOmOpEi = prOmOpE->at(i); - auto &prOeOpEi = prOeOpE->at(i); - for (int j = 0; j < 4; j++) { - auto prOeOpEij = prOmOpEi->at(j) + pAOmpE->at(j)->at(i)->timesFullColumn(rmem); - prOeOpEi->at(j) = prOeOpEij; - } - } - auto rpep = markerFrame->rpmp->plusFullColumn(markerFrame->aApm->timesFullColumn(rmem)); - pprOeOpEpE = EulerParameters::ppApEpEtimesColumn(rpep); - aAOe = aAOm->timesFullMatrix(aAme); - for (int i = 0; i < 4; i++) { - pAOepE->at(i) = pAOmpE->at(i)->timesFullMatrix(aAme); - } - auto aApe = markerFrame->aApm->timesFullMatrix(aAme); - ppAOepEpE = EulerParameters::ppApEpEtimesMatrix(aApe); - } - - void EndFrameqct::prePosIC() { - time = this->root()->mbdTimeValue(); - this->evalrmem(); - this->evalAme(); - EndFrameqc::prePosIC(); - } - - void EndFrameqct::evalrmem() { - if (rmemBlks) { - for (int i = 0; i < 3; i++) { - auto &expression = rmemBlks->at(i); - double value = expression->getValue(); - rmem->at(i) = value; - } - } - } - - void EndFrameqct::evalAme() { - if (phiThePsiBlks) { - auto phiThePsi = CREATE>::With(); - for (int i = 0; i < 3; i++) { - auto &expression = phiThePsiBlks->at(i); - auto value = expression->getValue(); - phiThePsi->at(i) = value; - } - phiThePsi->calc(); - aAme = phiThePsi->aA; - } - } - - void EndFrameqct::preVelIC() { - time = this->root()->mbdTimeValue(); - this->evalrmem(); - this->evalAme(); - Item::preVelIC(); - this->evalprmempt(); - this->evalpAmept(); - auto &aAOm = markerFrame->aAOm; - prOeOpt = aAOm->timesFullColumn(prmempt); - pAOept = aAOm->timesFullMatrix(pAmept); - } - - void EndFrameqct::postVelIC() { - auto &pAOmpE = markerFrame->pAOmpE; - for (int i = 0; i < 3; i++) { - auto &pprOeOpEpti = pprOeOpEpt->at(i); - for (int j = 0; j < 4; j++) { - auto pprOeOpEptij = pAOmpE->at(j)->at(i)->dot(prmempt); - pprOeOpEpti->atiput(j, pprOeOpEptij); - } - } - for (int i = 0; i < 4; i++) { - ppAOepEpt->atiput(i, pAOmpE->at(i)->timesFullMatrix(pAmept)); - } - } - - FColDsptr EndFrameqct::pAjOept(int j) { - return pAOept->column(j); - } - - FMatDsptr EndFrameqct::ppAjOepETpt(int jj) { - auto answer = std::make_shared(4, 3); - for (int i = 0; i < 4; i++) { - auto &answeri = answer->at(i); - auto &ppAOepEipt = ppAOepEpt->at(i); - for (int j = 0; j < 3; j++) { - auto &answerij = ppAOepEipt->at(j)->at(jj); - answeri->atiput(j, answerij); - } - } - return answer; - } - - FColDsptr EndFrameqct::ppAjOeptpt(int j) { - return ppAOeptpt->column(j); - } - - double EndFrameqct::priOeOpt(int i) { - return prOeOpt->at(i); - } - - FRowDsptr EndFrameqct::ppriOeOpEpt(int i) { - return pprOeOpEpt->at(i); - } - - double EndFrameqct::ppriOeOptpt(int i) { - return pprOeOptpt->at(i); - } - - void EndFrameqct::evalprmempt() { - if (rmemBlks) { - for (int i = 0; i < 3; i++) { - auto &derivative = prmemptBlks->at(i); - auto value = derivative->getValue(); - prmempt->at(i) = value; - } - } - } - - void EndFrameqct::evalpAmept() { - if (phiThePsiBlks) { - auto phiThePsi = CREATE>::With(); - auto phiThePsiDot = CREATE>::With(); - phiThePsiDot->phiThePsi = phiThePsi; - for (int i = 0; i < 3; i++) { - auto &expression = phiThePsiBlks->at(i); - auto &derivative = pPhiThePsiptBlks->at(i); - auto value = expression->getValue(); - auto valueDot = derivative->getValue(); - phiThePsi->at(i) = value; - phiThePsiDot->at(i) = valueDot; - } - phiThePsi->calc(); - phiThePsiDot->calc(); - pAmept = phiThePsiDot->aAdot; - } - } - - void EndFrameqct::evalpprmemptpt() { - if (rmemBlks) { - for (int i = 0; i < 3; i++) { - auto &secondDerivative = pprmemptptBlks->at(i); - auto value = secondDerivative->getValue(); - pprmemptpt->atiput(i, value); - } - } - } - - void EndFrameqct::evalppAmeptpt() { - if (phiThePsiBlks) { - auto phiThePsi = CREATE>::With(); - auto phiThePsiDot = CREATE>::With(); - phiThePsiDot->phiThePsi = phiThePsi; - auto phiThePsiDDot = CREATE>::With(); - phiThePsiDDot->phiThePsiDot = phiThePsiDot; - for (int i = 0; i < 3; i++) { - auto &expression = phiThePsiBlks->at(i); - auto &derivative = pPhiThePsiptBlks->at(i); - auto &secondDerivative = ppPhiThePsiptptBlks->at(i); - auto value = expression->getValue(); - auto valueDot = derivative->getValue(); - auto valueDDot = secondDerivative->getValue(); - phiThePsi->atiput(i, value); - phiThePsiDot->atiput(i, valueDot); - phiThePsiDDot->atiput(i, valueDDot); - } - phiThePsi->calc(); - phiThePsiDot->calc(); - phiThePsiDDot->calc(); - ppAmeptpt = phiThePsiDDot->aAddot; - } - } - - FColDsptr EndFrameqct::rmeO() { - return markerFrame->aAOm->timesFullColumn(rmem); - } - - FColDsptr EndFrameqct::rpep() { - auto &rpmp = markerFrame->rpmp; - auto &aApm = markerFrame->aApm; - auto rpep = rpmp->plusFullColumn(aApm->timesFullColumn(rmem)); - return rpep; - } - - void EndFrameqct::preAccIC() { - time = this->root()->mbdTimeValue(); - this->evalrmem(); - this->evalAme(); - Item::preVelIC(); - this->evalprmempt(); - this->evalpAmept(); - auto &aAOm = markerFrame->aAOm; - prOeOpt = aAOm->timesFullColumn(prmempt); - pAOept = aAOm->timesFullMatrix(pAmept); - Item::preAccIC(); - this->evalpprmemptpt(); - this->evalppAmeptpt(); - aAOm = markerFrame->aAOm; - pprOeOptpt = aAOm->timesFullColumn(pprmemptpt); - ppAOeptpt = aAOm->timesFullMatrix(ppAmeptpt); - } - bool EndFrameqct::isEndFrameqc() - { - return false; - } +EndFrameqct::EndFrameqct() { +} + +EndFrameqct::EndFrameqct(const char* str) : EndFrameqc(str) { +} + +void EndFrameqct::initialize() +{ + EndFrameqc::initialize(); + rmem = std::make_shared>(3); + prmempt = std::make_shared>(3); + pprmemptpt = std::make_shared>(3); + aAme = FullMatrix::identitysptr(3); + pAmept = std::make_shared>(3, 3); + ppAmeptpt = std::make_shared>(3, 3); + pprOeOpEpt = std::make_shared>(3, 4); + pprOeOptpt = std::make_shared>(3); + ppAOepEpt = std::make_shared>(4); + ppAOeptpt = std::make_shared>(3, 3); +} + +void EndFrameqct::initializeLocally() +{ + if (!rmemBlks) { + rmem->zeroSelf(); + prmempt->zeroSelf(); + pprmemptpt->zeroSelf(); + } + if (!phiThePsiBlks) { + aAme->identity(); + pAmept->zeroSelf(); + ppAmeptpt->zeroSelf(); + } +} + +void EndFrameqct::initializeGlobally() +{ + if (rmemBlks) { + initprmemptBlks(); + initpprmemptptBlks(); + } + if (phiThePsiBlks) { + initpPhiThePsiptBlks(); + initppPhiThePsiptptBlks(); + } +} + +void EndFrameqct::initprmemptBlks() +{ + auto& mbdTime = this->root()->time; + prmemptBlks = std::make_shared< FullColumn>(3); + for (int i = 0; i < 3; i++) { + auto& disp = rmemBlks->at(i); + auto var = disp->differentiateWRT(mbdTime); + auto vel = var->simplified(var); + prmemptBlks->at(i) = vel; + } +} + +void EndFrameqct::initpprmemptptBlks() +{ + auto& mbdTime = this->root()->time; + pprmemptptBlks = std::make_shared< FullColumn>(3); + for (int i = 0; i < 3; i++) { + auto& vel = prmemptBlks->at(i); + auto var = vel->differentiateWRT(mbdTime); + auto acc = var->simplified(var); + pprmemptptBlks->at(i) = acc; + } +} + +void EndFrameqct::initpPhiThePsiptBlks() +{ + auto& mbdTime = this->root()->time; + pPhiThePsiptBlks = std::make_shared< FullColumn>(3); + for (int i = 0; i < 3; i++) { + auto& angle = phiThePsiBlks->at(i); + auto var = angle->differentiateWRT(mbdTime); + //std::cout << "var " << *var << std::endl; + auto vel = var->simplified(var); + //std::cout << "vel " << *vel << std::endl; + pPhiThePsiptBlks->at(i) = vel; + } +} + +void EndFrameqct::initppPhiThePsiptptBlks() +{ + auto& mbdTime = this->root()->time; + ppPhiThePsiptptBlks = std::make_shared< FullColumn>(3); + for (int i = 0; i < 3; i++) { + auto& angleVel = pPhiThePsiptBlks->at(i); + auto var = angleVel->differentiateWRT(mbdTime); + auto angleAcc = var->simplified(var); + ppPhiThePsiptptBlks->at(i) = angleAcc; + } +} + +void EndFrameqct::postInput() +{ + this->evalrmem(); + this->evalAme(); + Item::postInput(); +} + +void EndFrameqct::calcPostDynCorrectorIteration() +{ + auto& rOmO = markerFrame->rOmO; + auto& aAOm = markerFrame->aAOm; + rOeO = rOmO->plusFullColumn(aAOm->timesFullColumn(rmem)); + auto& prOmOpE = markerFrame->prOmOpE; + auto& pAOmpE = markerFrame->pAOmpE; + for (int i = 0; i < 3; i++) + { + auto& prOmOpEi = prOmOpE->at(i); + auto& prOeOpEi = prOeOpE->at(i); + for (int j = 0; j < 4; j++) + { + auto prOeOpEij = prOmOpEi->at(j) + pAOmpE->at(j)->at(i)->timesFullColumn(rmem); + prOeOpEi->at(j) = prOeOpEij; + } + } + auto rpep = markerFrame->rpmp->plusFullColumn(markerFrame->aApm->timesFullColumn(rmem)); + pprOeOpEpE = EulerParameters::ppApEpEtimesColumn(rpep); + aAOe = aAOm->timesFullMatrix(aAme); + for (int i = 0; i < 4; i++) + { + pAOepE->at(i) = pAOmpE->at(i)->timesFullMatrix(aAme); + } + auto aApe = markerFrame->aApm->timesFullMatrix(aAme); + ppAOepEpE = EulerParameters::ppApEpEtimesMatrix(aApe); +} + +void EndFrameqct::prePosIC() +{ + time = this->root()->mbdTimeValue(); + this->evalrmem(); + this->evalAme(); + EndFrameqc::prePosIC(); +} + +void EndFrameqct::evalrmem() +{ + if (rmemBlks) { + for (int i = 0; i < 3; i++) + { + auto& expression = rmemBlks->at(i); + double value = expression->getValue(); + rmem->at(i) = value; + } + } +} + +void EndFrameqct::evalAme() +{ + if (phiThePsiBlks) { + auto phiThePsi = CREATE>::With(); + for (int i = 0; i < 3; i++) + { + auto& expression = phiThePsiBlks->at(i); + auto value = expression->getValue(); + phiThePsi->at(i) = value; + } + phiThePsi->calc(); + aAme = phiThePsi->aA; + } +} + +void EndFrameqct::preVelIC() +{ + time = this->root()->mbdTimeValue(); + this->evalrmem(); + this->evalAme(); + Item::preVelIC(); + this->evalprmempt(); + this->evalpAmept(); + auto& aAOm = markerFrame->aAOm; + prOeOpt = aAOm->timesFullColumn(prmempt); + pAOept = aAOm->timesFullMatrix(pAmept); +} + +void EndFrameqct::postVelIC() +{ + auto& pAOmpE = markerFrame->pAOmpE; + for (int i = 0; i < 3; i++) + { + auto& pprOeOpEpti = pprOeOpEpt->at(i); + for (int j = 0; j < 4; j++) + { + auto pprOeOpEptij = pAOmpE->at(j)->at(i)->dot(prmempt); + pprOeOpEpti->atiput(j, pprOeOpEptij); + } + } + for (int i = 0; i < 4; i++) + { + ppAOepEpt->atiput(i, pAOmpE->at(i)->timesFullMatrix(pAmept)); + } +} + +FColDsptr EndFrameqct::pAjOept(int j) +{ + return pAOept->column(j); +} + +FMatDsptr EndFrameqct::ppAjOepETpt(int jj) +{ + auto answer = std::make_shared>(4, 3); + for (int i = 0; i < 4; i++) + { + auto& answeri = answer->at(i); + auto& ppAOepEipt = ppAOepEpt->at(i); + for (int j = 0; j < 3; j++) + { + auto& answerij = ppAOepEipt->at(j)->at(jj); + answeri->atiput(j, answerij); + } + } + return answer; +} + +FColDsptr EndFrameqct::ppAjOeptpt(int j) +{ + return ppAOeptpt->column(j); +} + +double EndFrameqct::priOeOpt(int i) +{ + return prOeOpt->at(i); +} + +FRowDsptr EndFrameqct::ppriOeOpEpt(int i) +{ + return pprOeOpEpt->at(i); +} + +double EndFrameqct::ppriOeOptpt(int i) +{ + return pprOeOptpt->at(i); +} + +void EndFrameqct::evalprmempt() +{ + if (rmemBlks) { + for (int i = 0; i < 3; i++) + { + auto& derivative = prmemptBlks->at(i); + auto value = derivative->getValue(); + prmempt->at(i) = value; + } + } +} + +void EndFrameqct::evalpAmept() +{ + if (phiThePsiBlks) { + auto phiThePsi = CREATE>::With(); + auto phiThePsiDot = CREATE>::With(); + phiThePsiDot->phiThePsi = phiThePsi; + for (int i = 0; i < 3; i++) + { + auto& expression = phiThePsiBlks->at(i); + auto& derivative = pPhiThePsiptBlks->at(i); + auto value = expression->getValue(); + auto valueDot = derivative->getValue(); + phiThePsi->at(i) = value; + phiThePsiDot->at(i) = valueDot; + } + phiThePsi->calc(); + phiThePsiDot->calc(); + pAmept = phiThePsiDot->aAdot; + } +} + +void EndFrameqct::evalpprmemptpt() +{ + if (rmemBlks) { + for (int i = 0; i < 3; i++) + { + auto& secondDerivative = pprmemptptBlks->at(i); + auto value = secondDerivative->getValue(); + pprmemptpt->atiput(i, value); + } + } +} + +void EndFrameqct::evalppAmeptpt() +{ + if (phiThePsiBlks) { + auto phiThePsi = CREATE>::With(); + auto phiThePsiDot = CREATE>::With(); + phiThePsiDot->phiThePsi = phiThePsi; + auto phiThePsiDDot = CREATE>::With(); + phiThePsiDDot->phiThePsiDot = phiThePsiDot; + for (int i = 0; i < 3; i++) + { + auto& expression = phiThePsiBlks->at(i); + auto& derivative = pPhiThePsiptBlks->at(i); + auto& secondDerivative = ppPhiThePsiptptBlks->at(i); + auto value = expression->getValue(); + auto valueDot = derivative->getValue(); + auto valueDDot = secondDerivative->getValue(); + phiThePsi->atiput(i, value); + phiThePsiDot->atiput(i, valueDot); + phiThePsiDDot->atiput(i, valueDDot); + } + phiThePsi->calc(); + phiThePsiDot->calc(); + phiThePsiDDot->calc(); + ppAmeptpt = phiThePsiDDot->aAddot; + } +} + +FColDsptr EndFrameqct::rmeO() +{ + return markerFrame->aAOm->timesFullColumn(rmem); +} + +FColDsptr EndFrameqct::rpep() +{ + auto& rpmp = markerFrame->rpmp; + auto& aApm = markerFrame->aApm; + auto rpep = rpmp->plusFullColumn(aApm->timesFullColumn(rmem)); + return rpep; +} + +void EndFrameqct::preAccIC() +{ + time = this->root()->mbdTimeValue(); + this->evalrmem(); + this->evalAme(); + Item::preVelIC(); + this->evalprmempt(); + this->evalpAmept(); + auto& aAOm = markerFrame->aAOm; + prOeOpt = aAOm->timesFullColumn(prmempt); + pAOept = aAOm->timesFullMatrix(pAmept); + Item::preAccIC(); + this->evalpprmemptpt(); + this->evalppAmeptpt(); + aAOm = markerFrame->aAOm; + pprOeOptpt = aAOm->timesFullColumn(pprmemptpt); + ppAOeptpt = aAOm->timesFullMatrix(ppAmeptpt); +} + +bool MbD::EndFrameqct::isEndFrameqc() +{ + return false; } diff --git a/OndselSolver/EndFrameqct2.cpp b/OndselSolver/EndFrameqct2.cpp index b2f0ce4..3f9a86b 100644 --- a/OndselSolver/EndFrameqct2.cpp +++ b/OndselSolver/EndFrameqct2.cpp @@ -29,7 +29,8 @@ void EndFrameqct2::initpPhiThePsiptBlks() { auto& mbdTime = this->root()->time; auto eulerAngles = std::static_pointer_cast>(phiThePsiBlks); - pPhiThePsiptBlks = differentiateWRT(*eulerAngles, mbdTime); + //pPhiThePsiptBlks = differentiateWRT(*eulerAngles, mbdTime); + pPhiThePsiptBlks = eulerAngles->differentiateWRT(mbdTime); } void EndFrameqct2::initppPhiThePsiptptBlks() diff --git a/OndselSolver/EulerAngles.cpp b/OndselSolver/EulerAngles.cpp index 4c843e4..e2a88f9 100644 --- a/OndselSolver/EulerAngles.cpp +++ b/OndselSolver/EulerAngles.cpp @@ -7,76 +7,3 @@ ***************************************************************************/ #include "EulerAngles.h" - -namespace MbD { - template - inline void EulerAngles::initialize() - { - rotOrder = std::make_shared>(3); - rotOrder->at(0) = 0; - rotOrder->at(1) = 1; - rotOrder->at(2) = 2; - } - template<> - inline void EulerAngles::calc() - { - cA = std::make_shared>(3); - for (int i = 0; i < 3; i++) - { - auto axis = rotOrder->at(i); - auto angle = this->at(i)->getValue(); - if (axis == 1) { - cA->atiput(i, FullMatrixDouble::rotatex(angle)); - } - else if (axis == 2) { - cA->atiput(i, FullMatrixDouble::rotatey(angle)); - } - else if (axis == 3) { - cA->atiput(i, FullMatrixDouble::rotatez(angle)); - } - else { - throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats."); - } - } - aA = cA->at(0)->timesFullMatrix(cA->at(1)->timesFullMatrix(cA->at(2))); - } - template<> - inline void EulerAngles::calc() - { - cA = std::make_shared>(3); - for (int i = 0; i < 3; i++) - { - auto axis = rotOrder->at(i); - auto angle = this->at(i); - if (axis == 1) { - cA->atiput(i, FullMatrixDouble::rotatex(angle)); - } - else if (axis == 2) { - cA->atiput(i, FullMatrixDouble::rotatey(angle)); - } - else if (axis == 3) { - cA->atiput(i, FullMatrixDouble::rotatez(angle)); - } - else { - throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats."); - } - } - aA = cA->at(0)->timesFullMatrix(cA->at(1)->timesFullMatrix(cA->at(2))); - } - template - inline void EulerAngles::calc() - { - assert(false); - } - // type-specific helper functions - std::shared_ptr>> differentiateWRT(EulerAngles>& ref, std::shared_ptr var) - { - auto derivatives = std::make_shared>>(); - std::transform(ref.begin(), ref.end(), derivatives->begin(), - [var](std::shared_ptr term) { return term->differentiateWRT(var); } - ); - derivatives->aEulerAngles = &ref; - return derivatives; - } -} - diff --git a/OndselSolver/EulerAngles.h b/OndselSolver/EulerAngles.h index 216e6ad..a5e1502 100644 --- a/OndselSolver/EulerAngles.h +++ b/OndselSolver/EulerAngles.h @@ -16,6 +16,8 @@ #include "Symbolic.h" namespace MbD { + //template + //class EulerAnglesDot; template class EulerAngles : public EulerArray @@ -27,23 +29,87 @@ namespace MbD { EulerAngles(std::initializer_list list) : EulerArray{ list } {} void initialize() override; void calc() override; + std::shared_ptr> differentiateWRT(T var); void setRotOrder(int i, int j, int k); std::shared_ptr> rotOrder; FColFMatDsptr cA; FMatDsptr aA; - }; - template class EulerAngles>; - template class EulerAngles; - std::shared_ptr>> differentiateWRT(EulerAngles>& ref, std::shared_ptr var); - template - void EulerAngles::setRotOrder(int i, int j, int k) - { - rotOrder = std::make_shared>(3); - rotOrder->at(0) = i; - rotOrder->at(1) = j; - rotOrder->at(2) = k; - } + }; + template + inline void EulerAngles::initialize() + { + assert(false); + } + template<> + inline void EulerAngles::calc() + { + cA = std::make_shared>(3); + for (int i = 0; i < 3; i++) + { + auto axis = rotOrder->at(i); + auto angle = this->at(i)->getValue(); + if (axis == 1) { + cA->atiput(i, FullMatrix::rotatex(angle)); + } + else if (axis == 2) { + cA->atiput(i, FullMatrix::rotatey(angle)); + } + else if (axis == 3) { + cA->atiput(i, FullMatrix::rotatez(angle)); + } + else { + throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats."); + } + } + aA = cA->at(0)->timesFullMatrix(cA->at(1)->timesFullMatrix(cA->at(2))); + } + template<> + inline void EulerAngles::calc() + { + cA = std::make_shared>(3); + for (int i = 0; i < 3; i++) + { + auto axis = rotOrder->at(i); + auto angle = this->at(i); + if (axis == 1) { + cA->atiput(i, FullMatrix::rotatex(angle)); + } + else if (axis == 2) { + cA->atiput(i, FullMatrix::rotatey(angle)); + } + else if (axis == 3) { + cA->atiput(i, FullMatrix::rotatez(angle)); + } + else { + throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats."); + } + } + aA = cA->at(0)->timesFullMatrix(cA->at(1)->timesFullMatrix(cA->at(2))); + } + template + inline void EulerAngles::calc() + { + assert(false); + } + template + inline std::shared_ptr> EulerAngles::differentiateWRT(T var) + { + auto derivatives = std::make_shared>(); + std::transform(this->begin(), this->end(), derivatives->begin(), + [var](T term) { return term->differentiateWRT(var); } + ); + derivatives->aEulerAngles = this; + return derivatives; + } + template + inline void EulerAngles::setRotOrder(int i, int j, int k) + { + rotOrder = std::make_shared>(3); + rotOrder->at(0) = i; + rotOrder->at(1) = j; + rotOrder->at(2) = k; + } } diff --git a/OndselSolver/EulerAnglesDDot.h b/OndselSolver/EulerAnglesDDot.h index 44edb61..43f47be 100644 --- a/OndselSolver/EulerAnglesDDot.h +++ b/OndselSolver/EulerAnglesDDot.h @@ -47,13 +47,13 @@ namespace MbD { auto angleDot = aEulerAnglesDot->at(i)->getValue(); auto angleDDot = this->at(i)->getValue(); if (axis == 1) { - cAddot->atiput(i, FullMatrixDouble::rotatexrotDotrotDDot(angle, angleDot, angleDDot)); + cAddot->atiput(i, FullMatrix::rotatexrotDotrotDDot(angle, angleDot, angleDDot)); } else if (axis == 2) { - cAddot->atiput(i, FullMatrixDouble::rotateyrotDotrotDDot(angle, angleDot, angleDDot)); + cAddot->atiput(i, FullMatrix::rotateyrotDotrotDDot(angle, angleDot, angleDDot)); } else if (axis == 3) { - cAddot->atiput(i, FullMatrixDouble::rotatezrotDotrotDDot(angle, angleDot, angleDDot)); + cAddot->atiput(i, FullMatrix::rotatezrotDotrotDDot(angle, angleDot, angleDDot)); } else { throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats."); diff --git a/OndselSolver/EulerAnglesDot.h b/OndselSolver/EulerAnglesDot.h index baf755f..23a18cd 100644 --- a/OndselSolver/EulerAnglesDot.h +++ b/OndselSolver/EulerAnglesDot.h @@ -56,13 +56,13 @@ namespace MbD { auto angle = aEulerAngles->at(i)->getValue(); auto angleDot = this->at(i)->getValue(); if (axis == 1) { - cAdot->atiput(i, FullMatrixDouble::rotatexrotDot(angle, angleDot)); + cAdot->atiput(i, FullMatrix::rotatexrotDot(angle, angleDot)); } else if (axis == 2) { - cAdot->atiput(i, FullMatrixDouble::rotateyrotDot(angle, angleDot)); + cAdot->atiput(i, FullMatrix::rotateyrotDot(angle, angleDot)); } else if (axis == 3) { - cAdot->atiput(i, FullMatrixDouble::rotatezrotDot(angle, angleDot)); + cAdot->atiput(i, FullMatrix::rotatezrotDot(angle, angleDot)); } else { throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats."); diff --git a/OndselSolver/EulerAngleszxz.h b/OndselSolver/EulerAngleszxz.h index 23a7617..55dc48f 100644 --- a/OndselSolver/EulerAngleszxz.h +++ b/OndselSolver/EulerAngleszxz.h @@ -29,9 +29,9 @@ namespace MbD { template inline void EulerAngleszxz::initialize() { - phiA = FullMatrixDouble::identitysptr(3); - theA = FullMatrixDouble::identitysptr(3); - psiA = FullMatrixDouble::identitysptr(3); + phiA = FullMatrix::identitysptr(3); + theA = FullMatrix::identitysptr(3); + psiA = FullMatrix::identitysptr(3); } template inline void EulerAngleszxz::calc() diff --git a/OndselSolver/EulerAngleszxzDDot.h b/OndselSolver/EulerAngleszxzDDot.h index 10879af..69bc344 100644 --- a/OndselSolver/EulerAngleszxzDDot.h +++ b/OndselSolver/EulerAngleszxzDDot.h @@ -8,7 +8,6 @@ #pragma once -#include "FullMatrix.h" #include "EulerArray.h" #include "EulerAngleszxzDot.h" @@ -24,17 +23,16 @@ namespace MbD { void calc() override; std::shared_ptr> phiThePsiDot; - FMatDsptr phiAddot, theAddot, psiAddot; - std::shared_ptr aAddot; + FMatDsptr phiAddot, theAddot, psiAddot, aAddot; }; template inline void EulerAngleszxzDDot::initialize() { - phiAddot = std::make_shared(3, 3); + phiAddot = std::make_shared>(3, 3); phiAddot->zeroSelf(); - theAddot = std::make_shared(3, 3); + theAddot = std::make_shared>(3, 3); theAddot->zeroSelf(); - psiAddot = std::make_shared(3, 3); + psiAddot = std::make_shared>(3, 3); psiAddot->zeroSelf(); } template @@ -92,7 +90,7 @@ namespace MbD { + *(phiAdot->timesFullMatrix(theA->timesFullMatrix(psiAdot))) + *(phiA->timesFullMatrix(theAdot->timesFullMatrix(psiAdot))) + *(phiA->timesFullMatrix(theA->timesFullMatrix(psiAddot))); - aAddot = std::make_shared(mat); + aAddot = std::make_shared>(mat); } } diff --git a/OndselSolver/EulerAngleszxzDot.h b/OndselSolver/EulerAngleszxzDot.h index 03f56b7..2037e70 100644 --- a/OndselSolver/EulerAngleszxzDot.h +++ b/OndselSolver/EulerAngleszxzDot.h @@ -28,11 +28,11 @@ namespace MbD { template inline void EulerAngleszxzDot::initialize() { - phiAdot = std::make_shared(3, 3); + phiAdot = std::make_shared>(3, 3); phiAdot->zeroSelf(); - theAdot = std::make_shared(3, 3); + theAdot = std::make_shared>(3, 3); theAdot->zeroSelf(); - psiAdot = std::make_shared(3, 3); + psiAdot = std::make_shared>(3, 3); psiAdot->zeroSelf(); } template diff --git a/OndselSolver/EulerParameters.cpp b/OndselSolver/EulerParameters.cpp index dc31b61..a9cf51f 100644 --- a/OndselSolver/EulerParameters.cpp +++ b/OndselSolver/EulerParameters.cpp @@ -7,308 +7,5 @@ ***************************************************************************/ #include "EulerParameters.h" -#include "FullColumn.h" -#include "FullRow.h" -#include "FullMatrix.h" -namespace MbD { - template<> - FMatFColDsptr EulerParameters::ppApEpEtimesColumn(FColDsptr col) - { - double a2c0 = 2 * col->at(0); - double a2c1 = 2 * col->at(1); - double a2c2 = 2 * col->at(2); - double m2c0 = 0 - a2c0; - double m2c1 = 0 - a2c1; - double m2c2 = 0 - a2c2; - auto col00 = std::make_shared>(ListD{ a2c0, m2c1, m2c2 }); - auto col01 = std::make_shared>(ListD{ a2c1, a2c0, 0 }); - auto col02 = std::make_shared>(ListD{ a2c2, 0, a2c0 }); - auto col03 = std::make_shared>(ListD{ 0, m2c2, a2c1 }); - auto col11 = std::make_shared>(ListD{ m2c0, a2c1, m2c2 }); - auto col12 = std::make_shared>(ListD{ 0, a2c2, a2c1 }); - auto col13 = std::make_shared>(ListD{ a2c2, 0, m2c0 }); - auto col22 = std::make_shared>(ListD{ m2c0, m2c1, a2c2 }); - auto col23 = std::make_shared>(ListD{ m2c1, a2c0, 0 }); - auto col33 = std::make_shared>(ListD{ a2c0, a2c1, a2c2 }); - auto answer = std::make_shared(4, 4); - auto& row0 = answer->at(0); - row0->at(0) = col00; - row0->at(1) = col01; - row0->at(2) = col02; - row0->at(3) = col03; - auto& row1 = answer->at(1); - row1->at(0) = col01; - row1->at(1) = col11; - row1->at(2) = col12; - row1->at(3) = col13; - auto& row2 = answer->at(2); - row2->at(0) = col02; - row2->at(1) = col12; - row2->at(2) = col22; - row2->at(3) = col23; - auto& row3 = answer->at(3); - row3->at(0) = col03; - row3->at(1) = col13; - row3->at(2) = col23; - row3->at(3) = col33; - return answer; - } - - template<> - FMatDsptr EulerParameters::pCpEtimesColumn(FColDsptr col) - { - //"col size = 4." - auto c0 = col->at(0); - auto c1 = col->at(1); - auto c2 = col->at(2); - auto mc0 = -c0; - auto mc1 = -c1; - auto mc2 = -c2; - auto mc3 = -col->at(3); - auto answer = std::make_shared(3, 4); - auto& row0 = answer->at(0); - auto& row1 = answer->at(1); - auto& row2 = answer->at(2); - row0->atiput(0, mc3); - row0->atiput(1, mc2); - row0->atiput(2, c1); - row0->atiput(3, c0); - row1->atiput(0, c2); - row1->atiput(1, mc3); - row1->atiput(2, mc0); - row1->atiput(3, c1); - row2->atiput(0, mc1); - row2->atiput(1, c0); - row2->atiput(2, mc3); - row2->atiput(3, c2); - return answer; - } - - template - FMatDsptr EulerParameters::pCTpEtimesColumn(FColDsptr col) - { - //"col size = 3." - auto c0 = col->at(0); - auto c1 = col->at(1); - auto c2 = col->at(2); - auto mc0 = -c0; - auto mc1 = -c1; - auto mc2 = -c2; - auto answer = std::make_shared(4, 4); - auto& row0 = answer->at(0); - auto& row1 = answer->at(1); - auto& row2 = answer->at(2); - auto& row3 = answer->at(3); - row0->atiput(0, 0.0); - row0->atiput(1, c2); - row0->atiput(2, mc1); - row0->atiput(3, c0); - row1->atiput(0, mc2); - row1->atiput(1, 0.0); - row1->atiput(2, c0); - row1->atiput(3, c1); - row2->atiput(0, c1); - row2->atiput(1, mc0); - row2->atiput(2, 0.0); - row2->atiput(3, c2); - row3->atiput(0, mc0); - row3->atiput(1, mc1); - row3->atiput(2, mc2); - row3->atiput(3, 0.0); - return answer; - } - - template<> - FMatFMatDsptr EulerParameters::ppApEpEtimesMatrix(FMatDsptr mat) - { - FRowDsptr a2m0 = mat->at(0)->times(2.0); - FRowDsptr a2m1 = mat->at(1)->times(2.0); - FRowDsptr a2m2 = mat->at(2)->times(2.0); - FRowDsptr m2m0 = a2m0->negated(); - FRowDsptr m2m1 = a2m1->negated(); - FRowDsptr m2m2 = a2m2->negated(); - FRowDsptr zero = std::make_shared>(3, 0.0); - auto mat00 = std::make_shared(ListFRD{ a2m0, m2m1, m2m2 }); - auto mat01 = std::make_shared(ListFRD{ a2m1, a2m0, zero }); - auto mat02 = std::make_shared(ListFRD{ a2m2, zero, a2m0 }); - auto mat03 = std::make_shared(ListFRD{ zero, m2m2, a2m1 }); - auto mat11 = std::make_shared(ListFRD{ m2m0, a2m1, m2m2 }); - auto mat12 = std::make_shared(ListFRD{ zero, a2m2, a2m1 }); - auto mat13 = std::make_shared(ListFRD{ a2m2, zero, m2m0 }); - auto mat22 = std::make_shared(ListFRD{ m2m0, m2m1, a2m2 }); - auto mat23 = std::make_shared(ListFRD{ m2m1, a2m0, zero }); - auto mat33 = std::make_shared(ListFRD{ a2m0, a2m1, a2m2 }); - auto answer = std::make_shared(4, 4); - auto& row0 = answer->at(0); - row0->at(0) = mat00; - row0->at(1) = mat01; - row0->at(2) = mat02; - row0->at(3) = mat03; - auto& row1 = answer->at(1); - row1->at(0) = mat01; - row1->at(1) = mat11; - row1->at(2) = mat12; - row1->at(3) = mat13; - auto& row2 = answer->at(2); - row2->at(0) = mat02; - row2->at(1) = mat12; - row2->at(2) = mat22; - row2->at(3) = mat23; - auto& row3 = answer->at(3); - row3->at(0) = mat03; - row3->at(1) = mat13; - row3->at(2) = mat23; - row3->at(3) = mat33; - return answer; - } - - template // this is ALWAYS double; see note below. - void EulerParameters::initialize() - { - aA = FullMatrixDouble::identitysptr(3); - aB = std::make_shared(3, 4); - aC = std::make_shared(3, 4); - pApE = std::make_shared>(4); - for (int i = 0; i < 4; i++) - { - pApE->at(i) = std::make_shared(3, 3); - } - } - -// the following can't be valid as FullMatrix instatiatiates , yet -// this class needs to see the member functions of FullMatrix -//template<> -//void EulerParameters::initialize() -//{ -//} - - template - void EulerParameters::calc() - { - this->calcABC(); - this->calcpApE(); - } - template<> - void EulerParameters::calcABC() - { - double aE0 = this->at(0); - double aE1 = this->at(1); - double aE2 = this->at(2); - double aE3 = this->at(3); - double mE0 = -aE0; - double mE1 = -aE1; - double mE2 = -aE2; - FRowDsptr aBi; - aBi = aB->at(0); - aBi->at(0) = aE3; - aBi->at(1) = mE2; - aBi->at(2) = aE1; - aBi->at(3) = mE0; - aBi = aB->at(1); - aBi->at(0) = aE2; - aBi->at(1) = aE3; - aBi->at(2) = mE0; - aBi->at(3) = mE1; - aBi = aB->at(2); - aBi->at(0) = mE1; - aBi->at(1) = aE0; - aBi->at(2) = aE3; - aBi->at(3) = mE2; - FRowDsptr aCi; - aCi = aC->at(0); - aCi->at(0) = aE3; - aCi->at(1) = aE2; - aCi->at(2) = mE1; - aCi->at(3) = mE0; - aCi = aC->at(1); - aCi->at(0) = mE2; - aCi->at(1) = aE3; - aCi->at(2) = aE0; - aCi->at(3) = mE1; - aCi = aC->at(2); - aCi->at(0) = aE1; - aCi->at(1) = mE0; - aCi->at(2) = aE3; - aCi->at(3) = mE2; - - aA = aB->timesTransposeFullMatrix(aC); - } - template<> - void EulerParameters::calcpApE() - { - double a2E0 = 2.0 * (this->at(0)); - double a2E1 = 2.0 * (this->at(1)); - double a2E2 = 2.0 * (this->at(2)); - double a2E3 = 2.0 * (this->at(3)); - double m2E0 = -a2E0; - double m2E1 = -a2E1; - double m2E2 = -a2E2; - double m2E3 = -a2E3; - FMatDsptr pApEk; - pApEk = pApE->at(0); - FRowDsptr pAipEk; - pAipEk = pApEk->at(0); - pAipEk->at(0) = a2E0; - pAipEk->at(1) = a2E1; - pAipEk->at(2) = a2E2; - pAipEk = pApEk->at(1); - pAipEk->at(0) = a2E1; - pAipEk->at(1) = m2E0; - pAipEk->at(2) = m2E3; - pAipEk = pApEk->at(2); - pAipEk->at(0) = a2E2; - pAipEk->at(1) = a2E3; - pAipEk->at(2) = m2E0; - // - pApEk = pApE->at(1); - pAipEk = pApEk->at(0); - pAipEk->at(0) = m2E1; - pAipEk->at(1) = a2E0; - pAipEk->at(2) = a2E3; - pAipEk = pApEk->at(1); - pAipEk->at(0) = a2E0; - pAipEk->at(1) = a2E1; - pAipEk->at(2) = a2E2; - pAipEk = pApEk->at(2); - pAipEk->at(0) = m2E3; - pAipEk->at(1) = a2E2; - pAipEk->at(2) = m2E1; - // - pApEk = pApE->at(2); - pAipEk = pApEk->at(0); - pAipEk->at(0) = m2E2; - pAipEk->at(1) = m2E3; - pAipEk->at(2) = a2E0; - pAipEk = pApEk->at(1); - pAipEk->at(0) = a2E3; - pAipEk->at(1) = m2E2; - pAipEk->at(2) = a2E1; - pAipEk = pApEk->at(2); - pAipEk->at(0) = a2E0; - pAipEk->at(1) = a2E1; - pAipEk->at(2) = a2E2; - // - pApEk = pApE->at(3); - pAipEk = pApEk->at(0); - pAipEk->at(0) = a2E3; - pAipEk->at(1) = m2E2; - pAipEk->at(2) = a2E1; - pAipEk = pApEk->at(1); - pAipEk->at(0) = a2E2; - pAipEk->at(1) = a2E3; - pAipEk->at(2) = m2E0; - pAipEk = pApEk->at(2); - pAipEk->at(0) = m2E1; - pAipEk->at(1) = a2E0; - pAipEk->at(2) = a2E3; - } - template - void EulerParameters::conditionSelf() - { - EulerArray::conditionSelf(); - this->normalizeSelf(); - } - - template class EulerParameters; -} +using namespace MbD; diff --git a/OndselSolver/EulerParameters.h b/OndselSolver/EulerParameters.h index 35d32ae..94ff274 100644 --- a/OndselSolver/EulerParameters.h +++ b/OndselSolver/EulerParameters.h @@ -8,12 +8,9 @@ #pragma once -#include "FullColumn.ref.h" -#include "FullMatrix.ref.h" -#include "EulerParameters.ref.h" #include "EulerArray.h" -// #include "FullColumn.h" -// #include "FullMatrix.h" +#include "FullColumn.h" +#include "FullMatrix.h" namespace MbD { @@ -42,10 +39,11 @@ namespace MbD { this->initialize(); this->calc(); } - static std::shared_ptr ppApEpEtimesColumn(FColDsptr col); + + static std::shared_ptr>> ppApEpEtimesColumn(FColDsptr col); static FMatDsptr pCpEtimesColumn(FColDsptr col); static FMatDsptr pCTpEtimesColumn(FColDsptr col); - static std::shared_ptr ppApEpEtimesMatrix(FMatDsptr mat); + static std::shared_ptr>> ppApEpEtimesMatrix(FMatDsptr mat); void initialize() override; @@ -59,4 +57,295 @@ namespace MbD { FMatDsptr aC; FColFMatDsptr pApE; }; + + template<> + inline FMatFColDsptr EulerParameters::ppApEpEtimesColumn(FColDsptr col) + { + double a2c0 = 2 * col->at(0); + double a2c1 = 2 * col->at(1); + double a2c2 = 2 * col->at(2); + double m2c0 = 0 - a2c0; + double m2c1 = 0 - a2c1; + double m2c2 = 0 - a2c2; + auto col00 = std::make_shared>(ListD{ a2c0, m2c1, m2c2 }); + auto col01 = std::make_shared>(ListD{ a2c1, a2c0, 0 }); + auto col02 = std::make_shared>(ListD{ a2c2, 0, a2c0 }); + auto col03 = std::make_shared>(ListD{ 0, m2c2, a2c1 }); + auto col11 = std::make_shared>(ListD{ m2c0, a2c1, m2c2 }); + auto col12 = std::make_shared>(ListD{ 0, a2c2, a2c1 }); + auto col13 = std::make_shared>(ListD{ a2c2, 0, m2c0 }); + auto col22 = std::make_shared>(ListD{ m2c0, m2c1, a2c2 }); + auto col23 = std::make_shared>(ListD{ m2c1, a2c0, 0 }); + auto col33 = std::make_shared>(ListD{ a2c0, a2c1, a2c2 }); + auto answer = std::make_shared>(4, 4); + auto& row0 = answer->at(0); + row0->at(0) = col00; + row0->at(1) = col01; + row0->at(2) = col02; + row0->at(3) = col03; + auto& row1 = answer->at(1); + row1->at(0) = col01; + row1->at(1) = col11; + row1->at(2) = col12; + row1->at(3) = col13; + auto& row2 = answer->at(2); + row2->at(0) = col02; + row2->at(1) = col12; + row2->at(2) = col22; + row2->at(3) = col23; + auto& row3 = answer->at(3); + row3->at(0) = col03; + row3->at(1) = col13; + row3->at(2) = col23; + row3->at(3) = col33; + return answer; + } + + template<> + inline FMatDsptr EulerParameters::pCpEtimesColumn(FColDsptr col) + { + //"col size = 4." + auto c0 = col->at(0); + auto c1 = col->at(1); + auto c2 = col->at(2); + auto mc0 = -c0; + auto mc1 = -c1; + auto mc2 = -c2; + auto mc3 = -col->at(3); + auto answer = std::make_shared>(3, 4); + auto& row0 = answer->at(0); + auto& row1 = answer->at(1); + auto& row2 = answer->at(2); + row0->atiput(0, mc3); + row0->atiput(1, mc2); + row0->atiput(2, c1); + row0->atiput(3, c0); + row1->atiput(0, c2); + row1->atiput(1, mc3); + row1->atiput(2, mc0); + row1->atiput(3, c1); + row2->atiput(0, mc1); + row2->atiput(1, c0); + row2->atiput(2, mc3); + row2->atiput(3, c2); + return answer; + } + + template + inline FMatDsptr EulerParameters::pCTpEtimesColumn(FColDsptr col) + { + //"col size = 3." + auto c0 = col->at(0); + auto c1 = col->at(1); + auto c2 = col->at(2); + auto mc0 = -c0; + auto mc1 = -c1; + auto mc2 = -c2; + auto answer = std::make_shared>(4, 4); + auto& row0 = answer->at(0); + auto& row1 = answer->at(1); + auto& row2 = answer->at(2); + auto& row3 = answer->at(3); + row0->atiput(0, 0.0); + row0->atiput(1, c2); + row0->atiput(2, mc1); + row0->atiput(3, c0); + row1->atiput(0, mc2); + row1->atiput(1, 0.0); + row1->atiput(2, c0); + row1->atiput(3, c1); + row2->atiput(0, c1); + row2->atiput(1, mc0); + row2->atiput(2, 0.0); + row2->atiput(3, c2); + row3->atiput(0, mc0); + row3->atiput(1, mc1); + row3->atiput(2, mc2); + row3->atiput(3, 0.0); + return answer; + } + + template<> + inline FMatFMatDsptr EulerParameters::ppApEpEtimesMatrix(FMatDsptr mat) + { + FRowDsptr a2m0 = mat->at(0)->times(2.0); + FRowDsptr a2m1 = mat->at(1)->times(2.0); + FRowDsptr a2m2 = mat->at(2)->times(2.0); + FRowDsptr m2m0 = a2m0->negated(); + FRowDsptr m2m1 = a2m1->negated(); + FRowDsptr m2m2 = a2m2->negated(); + FRowDsptr zero = std::make_shared>(3, 0.0); + auto mat00 = std::make_shared>(ListFRD{ a2m0, m2m1, m2m2 }); + auto mat01 = std::make_shared>(ListFRD{ a2m1, a2m0, zero }); + auto mat02 = std::make_shared>(ListFRD{ a2m2, zero, a2m0 }); + auto mat03 = std::make_shared>(ListFRD{ zero, m2m2, a2m1 }); + auto mat11 = std::make_shared>(ListFRD{ m2m0, a2m1, m2m2 }); + auto mat12 = std::make_shared>(ListFRD{ zero, a2m2, a2m1 }); + auto mat13 = std::make_shared>(ListFRD{ a2m2, zero, m2m0 }); + auto mat22 = std::make_shared>(ListFRD{ m2m0, m2m1, a2m2 }); + auto mat23 = std::make_shared>(ListFRD{ m2m1, a2m0, zero }); + auto mat33 = std::make_shared>(ListFRD{ a2m0, a2m1, a2m2 }); + auto answer = std::make_shared>(4, 4); + auto& row0 = answer->at(0); + row0->at(0) = mat00; + row0->at(1) = mat01; + row0->at(2) = mat02; + row0->at(3) = mat03; + auto& row1 = answer->at(1); + row1->at(0) = mat01; + row1->at(1) = mat11; + row1->at(2) = mat12; + row1->at(3) = mat13; + auto& row2 = answer->at(2); + row2->at(0) = mat02; + row2->at(1) = mat12; + row2->at(2) = mat22; + row2->at(3) = mat23; + auto& row3 = answer->at(3); + row3->at(0) = mat03; + row3->at(1) = mat13; + row3->at(2) = mat23; + row3->at(3) = mat33; + return answer; + } + + template<> + inline void EulerParameters::initialize() + { + aA = FullMatrix::identitysptr(3); + aB = std::make_shared>(3, 4); + aC = std::make_shared>(3, 4); + pApE = std::make_shared>(4); + for (int i = 0; i < 4; i++) + { + pApE->at(i) = std::make_shared>(3, 3); + } + } + template + inline void EulerParameters::calc() + { + this->calcABC(); + this->calcpApE(); + } + template<> + inline void EulerParameters::calcABC() + { + double aE0 = this->at(0); + double aE1 = this->at(1); + double aE2 = this->at(2); + double aE3 = this->at(3); + double mE0 = -aE0; + double mE1 = -aE1; + double mE2 = -aE2; + FRowDsptr aBi; + aBi = aB->at(0); + aBi->at(0) = aE3; + aBi->at(1) = mE2; + aBi->at(2) = aE1; + aBi->at(3) = mE0; + aBi = aB->at(1); + aBi->at(0) = aE2; + aBi->at(1) = aE3; + aBi->at(2) = mE0; + aBi->at(3) = mE1; + aBi = aB->at(2); + aBi->at(0) = mE1; + aBi->at(1) = aE0; + aBi->at(2) = aE3; + aBi->at(3) = mE2; + FRowDsptr aCi; + aCi = aC->at(0); + aCi->at(0) = aE3; + aCi->at(1) = aE2; + aCi->at(2) = mE1; + aCi->at(3) = mE0; + aCi = aC->at(1); + aCi->at(0) = mE2; + aCi->at(1) = aE3; + aCi->at(2) = aE0; + aCi->at(3) = mE1; + aCi = aC->at(2); + aCi->at(0) = aE1; + aCi->at(1) = mE0; + aCi->at(2) = aE3; + aCi->at(3) = mE2; + + aA = aB->timesTransposeFullMatrix(aC); + } + template<> + inline void EulerParameters::calcpApE() + { + double a2E0 = 2.0 * (this->at(0)); + double a2E1 = 2.0 * (this->at(1)); + double a2E2 = 2.0 * (this->at(2)); + double a2E3 = 2.0 * (this->at(3)); + double m2E0 = -a2E0; + double m2E1 = -a2E1; + double m2E2 = -a2E2; + double m2E3 = -a2E3; + FMatDsptr pApEk; + pApEk = pApE->at(0); + FRowDsptr pAipEk; + pAipEk = pApEk->at(0); + pAipEk->at(0) = a2E0; + pAipEk->at(1) = a2E1; + pAipEk->at(2) = a2E2; + pAipEk = pApEk->at(1); + pAipEk->at(0) = a2E1; + pAipEk->at(1) = m2E0; + pAipEk->at(2) = m2E3; + pAipEk = pApEk->at(2); + pAipEk->at(0) = a2E2; + pAipEk->at(1) = a2E3; + pAipEk->at(2) = m2E0; + // + pApEk = pApE->at(1); + pAipEk = pApEk->at(0); + pAipEk->at(0) = m2E1; + pAipEk->at(1) = a2E0; + pAipEk->at(2) = a2E3; + pAipEk = pApEk->at(1); + pAipEk->at(0) = a2E0; + pAipEk->at(1) = a2E1; + pAipEk->at(2) = a2E2; + pAipEk = pApEk->at(2); + pAipEk->at(0) = m2E3; + pAipEk->at(1) = a2E2; + pAipEk->at(2) = m2E1; + // + pApEk = pApE->at(2); + pAipEk = pApEk->at(0); + pAipEk->at(0) = m2E2; + pAipEk->at(1) = m2E3; + pAipEk->at(2) = a2E0; + pAipEk = pApEk->at(1); + pAipEk->at(0) = a2E3; + pAipEk->at(1) = m2E2; + pAipEk->at(2) = a2E1; + pAipEk = pApEk->at(2); + pAipEk->at(0) = a2E0; + pAipEk->at(1) = a2E1; + pAipEk->at(2) = a2E2; + // + pApEk = pApE->at(3); + pAipEk = pApEk->at(0); + pAipEk->at(0) = a2E3; + pAipEk->at(1) = m2E2; + pAipEk->at(2) = a2E1; + pAipEk = pApEk->at(1); + pAipEk->at(0) = a2E2; + pAipEk->at(1) = a2E3; + pAipEk->at(2) = m2E0; + pAipEk = pApEk->at(2); + pAipEk->at(0) = m2E1; + pAipEk->at(1) = a2E0; + pAipEk->at(2) = a2E3; + } + template + inline void EulerParameters::conditionSelf() + { + EulerArray::conditionSelf(); + this->normalizeSelf(); + } } + diff --git a/OndselSolver/EulerParametersDot.h b/OndselSolver/EulerParametersDot.h index 283408c..e4da335 100644 --- a/OndselSolver/EulerParametersDot.h +++ b/OndselSolver/EulerParametersDot.h @@ -55,13 +55,13 @@ namespace MbD { template inline void EulerParametersDot::initialize() { - aAdot = std::make_shared(3, 3); - aBdot = std::make_shared(3, 4); - aCdot = std::make_shared(3, 4); + aAdot = std::make_shared>(3, 3); + aBdot = std::make_shared>(3, 4); + aCdot = std::make_shared>(3, 4); pAdotpE = std::make_shared>(4); for (int i = 0; i < 4; i++) { - pAdotpE->at(i) = std::make_shared(3, 3); + pAdotpE->at(i) = std::make_shared>(3, 3); } } diff --git a/OndselSolver/FullColumn.cpp b/OndselSolver/FullColumn.cpp index d77458b..ae0ee14 100644 --- a/OndselSolver/FullColumn.cpp +++ b/OndselSolver/FullColumn.cpp @@ -12,193 +12,4 @@ #include "FullRow.h" #include "FullMatrix.h" -namespace MbD { - template - FColsptr FullColumn::plusFullColumn(FColsptr fullCol) - { - int n = (int)this->size(); - auto answer = std::make_shared>(n); - for (int i = 0; i < n; i++) { - answer->at(i) = this->at(i) + fullCol->at(i); - } - return answer; - } - template - FColsptr FullColumn::plusFullColumntimes(FColsptr fullCol, double factor) - { - assert(false); - return FColsptr(); - } - template - FColsptr FullColumn::minusFullColumn(FColsptr fullCol) - { - int n = (int)this->size(); - auto answer = std::make_shared>(n); - for (int i = 0; i < n; i++) { - answer->at(i) = this->at(i) - fullCol->at(i); - } - return answer; - } - template<> - FColDsptr FullColumn::times(double a) - { - int n = (int)this->size(); - auto answer = std::make_shared>(n); - for (int i = 0; i < n; i++) { - answer->at(i) = this->at(i) * a; - } - return answer; - } - template - FColsptr FullColumn::times(T a) - { - assert(false); - auto answer = std::make_shared>(); - return answer; - } - template - FColsptr FullColumn::negated() - { - return this->times((T) - 1.0); - } - template - void FullColumn::atiputFullColumn(int i, FColsptr fullCol) - { - for (int ii = 0; ii < fullCol->size(); ii++) - { - this->at(i + ii) = fullCol->at(ii); - } - } - template - void FullColumn::atiplusFullColumn(int i, FColsptr fullCol) - { - for (int ii = 0; ii < fullCol->size(); ii++) - { - this->at(i + ii) += fullCol->at(ii); - } - } - template - void FullColumn::equalSelfPlusFullColumnAt(FColsptr fullCol, int ii) - { - //self is subcolumn of fullCol - for (int i = 0; i < this->size(); i++) - { - this->at(i) += fullCol->at(ii + i); - } - } - template - void FullColumn::atiminusFullColumn(int i1, FColsptr fullCol) - { - for (int ii = 0; ii < fullCol->size(); ii++) - { - int i = i1 + ii; - this->at(i) -= fullCol->at(ii); - } - } - template - void FullColumn::equalFullColumnAt(FColsptr fullCol, int i) - { - this->equalArrayAt(fullCol, i); - //for (int ii = 0; ii < this->size(); ii++) - //{ - // this->at(ii) = fullCol->at(i + ii); - //} - } - template<> - FColDsptr FullColumn::copy() - { - auto n = (int)this->size(); - auto answer = std::make_shared>(n); - for (int i = 0; i < n; i++) - { - answer->at(i) = this->at(i); - } - return answer; - } - template - FRowsptr FullColumn::transpose() - { - return std::make_shared>(*this); - } - template - void FullColumn::atiplusFullColumntimes(int i1, FColsptr fullCol, T factor) - { - for (int ii = 0; ii < fullCol->size(); ii++) - { - int i = i1 + ii; - this->at(i) += fullCol->at(ii) * factor; - } - } - template - T FullColumn::transposeTimesFullColumn(const FColsptr fullCol) - { - return this->dotVec(fullCol); - } - template - void FullColumn::equalSelfPlusFullColumntimes(FColsptr fullCol, T factor) - { - this->equalSelfPlusFullVectortimes(fullCol, factor); - } - template - FColsptr FullColumn::cross(FColsptr fullCol) - { - auto a0 = this->at(0); - auto a1 = this->at(1); - auto a2 = this->at(2); - auto b0 = fullCol->at(0); - auto b1 = fullCol->at(1); - auto b2 = fullCol->at(2); - auto answer = std::make_shared>(3); - answer->atiput(0, a1 * b2 - (a2 * b1)); - answer->atiput(1, a2 * b0 - (a0 * b2)); - answer->atiput(2, a0 * b1 - (a1 * b0)); - return answer; - } - //template<> - //inline std::shared_ptr> FullColumn::simplified() - //{ - // auto n = (int)this->size(); - // auto answer = std::make_shared>(n); - // for (int i = 0; i < n; i++) - // { - // auto func = this->at(i); - // answer->at(i) = func->simplified(func); - // } - // return answer; - //} - template - FColsptr FullColumn::simplified() - { -// assert(false); - return FColsptr(); - } - // instantiate on purpose to make visible in library api: - template class FullColumn; - template class FullColumn; - template - double FullColumn::dot(std::shared_ptr> vec) - { - int n = (int)this->size(); - double answer = 0.0; - for (int i = 0; i < n; i++) { - answer += this->at(i) * vec->at(i); - } - return answer; - } - template - std::shared_ptr> FullColumn::dot(std::shared_ptr>>> vecvec) - { - int ncol = (int)this->size(); - auto nelem = vecvec->at(0)->size(); - auto answer = std::make_shared>(nelem); - for (int k = 0; k < nelem; k++) { - auto sum = (T)0; - for (int i = 0; i < ncol; i++) - { - sum += this->at(i) * vecvec->at(i)->at(k); - } - answer->at(k) = sum; - } - return answer; - } -} \ No newline at end of file +using namespace MbD; diff --git a/OndselSolver/FullColumn.h b/OndselSolver/FullColumn.h index 3bc2359..a3a06a7 100644 --- a/OndselSolver/FullColumn.h +++ b/OndselSolver/FullColumn.h @@ -5,19 +5,24 @@ * * * See LICENSE file for details about copyright. * ***************************************************************************/ - + #pragma once #include #include #include "FullVector.h" -#include "FullColumn.ref.h" -#include "FullRow.ref.h" -#include "FullMatrix.ref.h" -//#include "Symbolic.h" namespace MbD { + template + class FullColumn; + using FColDsptr = std::shared_ptr>; + template + using FColsptr = std::shared_ptr>; + template + class FullRow; + template + using FRowsptr = std::shared_ptr>; class Symbolic; template @@ -31,7 +36,6 @@ namespace MbD { FullColumn(typename std::vector::iterator begin, typename std::vector::iterator end) : FullVector(begin, end) {} FullColumn(std::initializer_list list) : FullVector{ list } {} FColsptr plusFullColumn(FColsptr fullCol); - FColsptr plusFullColumntimes(FColsptr fullCol, double factor); FColsptr minusFullColumn(FColsptr fullCol); FColsptr times(T a); FColsptr negated(); @@ -43,44 +47,203 @@ namespace MbD { FColsptr copy(); FRowsptr transpose(); void atiplusFullColumntimes(int i, FColsptr fullCol, T factor); - T transposeTimesFullColumn(const FColsptr fullCol); + T transposeTimesFullColumn(const FColsptr fullCol); void equalSelfPlusFullColumntimes(FColsptr fullCol, T factor); FColsptr cross(FColsptr fullCol); FColsptr simplified(); double dot(std::shared_ptr> vec); - std::shared_ptr> cloneFcSptr(); - double dotVec(std::shared_ptr> vec); std::shared_ptr> dot(std::shared_ptr>>> vecvec); std::ostream& printOn(std::ostream& s) const override; }; - // the following "printOn" needs to be in the header for unknown reasons linker - template - std::ostream& FullColumn::printOn(std::ostream& s) const - { - s << "FullCol{"; - s << this->at(0); - for (int i = 1; i < this->size(); i++) - { - s << ", " << this->at(i); - } - s << "}"; - return s; - } - template - std::shared_ptr> FullColumn::cloneFcSptr() - { - return std::make_shared>(*this); - } - template - double FullColumn::dotVec(std::shared_ptr> vec) - { - int n = (int)this->size(); - double answer = 0.0; - for (int i = 0; i < n; i++) { - answer += this->at(i) * vec->at(i); - } - return answer; - } - + template + inline FColsptr FullColumn::plusFullColumn(FColsptr fullCol) + { + int n = (int) this->size(); + auto answer = std::make_shared>(n); + for (int i = 0; i < n; i++) { + answer->at(i) = this->at(i) + fullCol->at(i); + } + return answer; + } + template + inline FColsptr FullColumn::minusFullColumn(FColsptr fullCol) + { + int n = (int) this->size(); + auto answer = std::make_shared>(n); + for (int i = 0; i < n; i++) { + answer->at(i) = this->at(i) - fullCol->at(i); + } + return answer; + } + template<> + inline FColDsptr FullColumn::times(double a) + { + int n = (int)this->size(); + auto answer = std::make_shared>(n); + for (int i = 0; i < n; i++) { + answer->at(i) = this->at(i) * a; + } + return answer; + } + template + inline FColsptr FullColumn::times(T a) + { + assert(false); + } + template + inline FColsptr FullColumn::negated() + { + return this->times(-1.0); + } + template + inline void FullColumn::atiputFullColumn(int i, FColsptr fullCol) + { + for (int ii = 0; ii < fullCol->size(); ii++) + { + this->at(i + ii) = fullCol->at(ii); + } + } + template + inline void FullColumn::atiplusFullColumn(int i, FColsptr fullCol) + { + for (int ii = 0; ii < fullCol->size(); ii++) + { + this->at(i + ii) += fullCol->at(ii); + } + } + template + inline void FullColumn::equalSelfPlusFullColumnAt(FColsptr fullCol, int ii) + { + //self is subcolumn of fullCol + for (int i = 0; i < this->size(); i++) + { + this->at(i) += fullCol->at(ii + i); + } + } + template + inline void FullColumn::atiminusFullColumn(int i1, FColsptr fullCol) + { + for (int ii = 0; ii < fullCol->size(); ii++) + { + int i = i1 + ii; + this->at(i) -= fullCol->at(ii); + } + } + template + inline void FullColumn::equalFullColumnAt(FColsptr fullCol, int i) + { + this->equalArrayAt(fullCol, i); + //for (int ii = 0; ii < this->size(); ii++) + //{ + // this->at(ii) = fullCol->at(i + ii); + //} + } + template<> + inline FColDsptr FullColumn::copy() + { + auto n = (int) this->size(); + auto answer = std::make_shared>(n); + for (int i = 0; i < n; i++) + { + answer->at(i) = this->at(i); + } + return answer; + } + template + inline FRowsptr FullColumn::transpose() + { + return std::make_shared>(*this); + } + template + inline void FullColumn::atiplusFullColumntimes(int i1, FColsptr fullCol, T factor) + { + for (int ii = 0; ii < fullCol->size(); ii++) + { + int i = i1 + ii; + this->at(i) += fullCol->at(ii) * factor; + } + } + template + inline T FullColumn::transposeTimesFullColumn(const FColsptr fullCol) + { + return this->dot(fullCol); + } + template + inline void FullColumn::equalSelfPlusFullColumntimes(FColsptr fullCol, T factor) + { + this->equalSelfPlusFullVectortimes(fullCol, factor); + } + template + inline FColsptr FullColumn::cross(FColsptr fullCol) + { + auto a0 = this->at(0); + auto a1 = this->at(1); + auto a2 = this->at(2); + auto b0 = fullCol->at(0); + auto b1 = fullCol->at(1); + auto b2 = fullCol->at(2); + auto answer = std::make_shared>(3); + answer->atiput(0, a1 * b2 - (a2 * b1)); + answer->atiput(1, a2 * b0 - (a0 * b2)); + answer->atiput(2, a0 * b1 - (a1 * b0)); + return answer; + } + //template<> + //inline std::shared_ptr> FullColumn::simplified() + //{ + // auto n = this->size(); + // auto answer = std::make_shared>(n); + // for (int i = 0; i < n; i++) + // { + // auto func = this->at(i); + // answer->at(i) = func->simplified(func); + // } + // return answer; + //} + template + inline FColsptr FullColumn::simplified() + { + assert(false); + return FColsptr(); + } + template + inline double FullColumn::dot(std::shared_ptr> vec) + { + int n = (int)this->size(); + double answer = 0.0; + for (int i = 0; i < n; i++) { + answer += this->at(i) * vec->at(i); + } + return answer; + } + template + inline std::shared_ptr> FullColumn::dot(std::shared_ptr>>> vecvec) + { + int ncol = (int)this->size(); + auto nelem = vecvec->at(0)->size(); + auto answer = std::make_shared>(nelem); + for (int k = 0; k < nelem; k++) { + auto sum = 0.0; + for (int i = 0; i < ncol; i++) + { + sum += this->at(i) * vecvec->at(i)->at(k); + } + answer->at(k) = sum; + } + return answer; + } + template + inline std::ostream& FullColumn::printOn(std::ostream& s) const + { + s << "FullCol{"; + s << this->at(0); + for (int i = 1; i < this->size(); i++) + { + s << ", " << this->at(i); + } + s << "}"; + return s; + } } + diff --git a/OndselSolver/FullColumn.ref.h b/OndselSolver/FullColumn.ref.h deleted file mode 100644 index d622179..0000000 --- a/OndselSolver/FullColumn.ref.h +++ /dev/null @@ -1,14 +0,0 @@ -#pragma once - -#include -#include - -namespace MbD { - template - class FullColumn; - - using FColDsptr = std::shared_ptr>; - - template - using FColsptr = std::shared_ptr>; -} diff --git a/OndselSolver/FullMatrix.cpp b/OndselSolver/FullMatrix.cpp index f338e70..65b2c5a 100644 --- a/OndselSolver/FullMatrix.cpp +++ b/OndselSolver/FullMatrix.cpp @@ -7,678 +7,6 @@ ***************************************************************************/ #include "FullMatrix.h" -#include "FullColumn.h" -#include "FullRow.h" -#include "DiagonalMatrix.h" -#include "EulerParameters.h" -namespace MbD { - FColsptr FullMatrixDouble::timesFullColumn(FColsptr fullCol) - { - return this->timesFullColumn(fullCol.get()); -// auto nrow = this->nrow(); -// auto answer = std::make_shared>(nrow); -// for (int i = 0; i < nrow; i++) -// { -// answer->at(i) = this->at(i)->timesFullColumn(fullCol); -// } -// return answer; - } - FColsptr FullMatrixDouble::timesFullColumn(FullColumn* fullCol) // local - { - //"a*b = a(i,j)b(j) sum j." - auto nrow = this->nrow(); - auto answer = std::make_shared>(nrow); - for (int i = 0; i < nrow; i++) - { - answer->at(i) = this->at(i)->timesFullColumn(fullCol); - } - return answer; - } - std::shared_ptr FullMatrixDouble::rotatex(double the) - { - auto sthe = std::sin(the); - auto cthe = std::cos(the); - auto rotMat = std::make_shared(3, 3); - auto row0 = rotMat->at(0); - row0->atiput(0, 1.0); - row0->atiput(1, 0.0); - row0->atiput(2, 0.0); - auto row1 = rotMat->at(1); - row1->atiput(0, 0.0); - row1->atiput(1, cthe); - row1->atiput(2, -sthe); - auto row2 = rotMat->at(2); - row2->atiput(0, 0.0); - row2->atiput(1, sthe); - row2->atiput(2, cthe); - return rotMat; - } - std::shared_ptr FullMatrixDouble::rotatey(double the) - { - auto sthe = std::sin(the); - auto cthe = std::cos(the); - auto rotMat = std::make_shared(3, 3); - auto row0 = rotMat->at(0); - row0->atiput(0, cthe); - row0->atiput(1, 0.0); - row0->atiput(2, sthe); - auto row1 = rotMat->at(1); - row1->atiput(0, 0.0); - row1->atiput(1, 1.0); - row1->atiput(2, 0.0); - auto row2 = rotMat->at(2); - row2->atiput(0, -sthe); - row2->atiput(1, 0.0); - row2->atiput(2, cthe); - return rotMat; - } - std::shared_ptr FullMatrixDouble::rotatez(double the) - { - auto sthe = std::sin(the); - auto cthe = std::cos(the); - auto rotMat = std::make_shared(3, 3); - auto row0 = rotMat->at(0); - row0->atiput(0, cthe); - row0->atiput(1, -sthe); - row0->atiput(2, 0.0); - auto row1 = rotMat->at(1); - row1->atiput(0, sthe); - row1->atiput(1, cthe); - row1->atiput(2, 0.0); - auto row2 = rotMat->at(2); - row2->atiput(0, 0.0); - row2->atiput(1, 0.0); - row2->atiput(2, 1.0); - return rotMat; - } - std::shared_ptr FullMatrixDouble::rotatexrotDot(double the, double thedot) - { - auto sthe = std::sin(the); - auto cthe = std::cos(the); - auto sthedot = cthe * thedot; - auto cthedot = -sthe * thedot; - auto rotMat = std::make_shared(3, 3); - auto row0 = rotMat->at(0); - row0->atiput(0, 0.0); - row0->atiput(1, 0.0); - row0->atiput(2, 0.0); - auto row1 = rotMat->at(1); - row1->atiput(0, 0.0); - row1->atiput(1, cthedot); - row1->atiput(2, -sthedot); - auto row2 = rotMat->at(2); - row2->atiput(0, 0.0); - row2->atiput(1, sthedot); - row2->atiput(2, cthedot); - return rotMat; - } - std::shared_ptr FullMatrixDouble::rotateyrotDot(double the, double thedot) - { - auto sthe = std::sin(the); - auto cthe = std::cos(the); - auto sthedot = cthe * thedot; - auto cthedot = -sthe * thedot; - auto rotMat = std::make_shared(3, 3); - auto row0 = rotMat->at(0); - row0->atiput(0, cthedot); - row0->atiput(1, 0.0); - row0->atiput(2, sthedot); - auto row1 = rotMat->at(1); - row1->atiput(0, 0.0); - row1->atiput(1, 0.0); - row1->atiput(2, 0.0); - auto row2 = rotMat->at(2); - row2->atiput(0, -sthedot); - row2->atiput(1, 0.0); - row2->atiput(2, cthedot); - return rotMat; - } - std::shared_ptr FullMatrixDouble::rotatezrotDot(double the, double thedot) - { - auto sthe = std::sin(the); - auto cthe = std::cos(the); - auto sthedot = cthe * thedot; - auto cthedot = -sthe * thedot; - auto rotMat = std::make_shared(3, 3); - auto row0 = rotMat->at(0); - row0->atiput(0, cthedot); - row0->atiput(1, -sthedot); - row0->atiput(2, 0.0); - auto row1 = rotMat->at(1); - row1->atiput(0, sthedot); - row1->atiput(1, cthedot); - row1->atiput(2, 0.0); - auto row2 = rotMat->at(2); - row2->atiput(0, 0.0); - row2->atiput(1, 0.0); - row2->atiput(2, 0.0); - return rotMat; - } - std::shared_ptr FullMatrixDouble::rotatexrotDotrotDDot(double the, double thedot, double theddot) - { - auto sthe = std::sin(the); - auto cthe = std::cos(the); - auto sthedot = cthe * thedot; - auto cthedot = -sthe * thedot; - auto stheddot = cthedot * thedot + (cthe * theddot); - auto ctheddot = -(sthedot * thedot) - (sthe * theddot); - auto rotMat = std::make_shared(3, 3); - auto row0 = rotMat->at(0); - row0->atiput(0, 0.0); - row0->atiput(1, 0.0); - row0->atiput(2, 0.0); - auto row1 = rotMat->at(1); - row1->atiput(0, 0.0); - row1->atiput(1, ctheddot); - row1->atiput(2, -stheddot); - auto row2 = rotMat->at(2); - row2->atiput(0, 0.0); - row2->atiput(1, stheddot); - row2->atiput(2, ctheddot); - return rotMat; - } - std::shared_ptr FullMatrixDouble::rotateyrotDotrotDDot(double the, double thedot, double theddot) - { - auto sthe = std::sin(the); - auto cthe = std::cos(the); - auto sthedot = cthe * thedot; - auto cthedot = -sthe * thedot; - auto stheddot = cthedot * thedot + (cthe * theddot); - auto ctheddot = -(sthedot * thedot) - (sthe * theddot); - auto rotMat = std::make_shared(3, 3); - auto row0 = rotMat->at(0); - row0->atiput(0, ctheddot); - row0->atiput(1, 0.0); - row0->atiput(2, stheddot); - auto row1 = rotMat->at(1); - row1->atiput(0, 0.0); - row1->atiput(1, 0.0); - row1->atiput(2, 0.0); - auto row2 = rotMat->at(2); - row2->atiput(0, -stheddot); - row2->atiput(1, 0.0); - row2->atiput(2, ctheddot); - return rotMat; - } - std::shared_ptr FullMatrixDouble::rotatezrotDotrotDDot(double the, double thedot, double theddot) - { - auto sthe = std::sin(the); - auto cthe = std::cos(the); - auto sthedot = cthe * thedot; - auto cthedot = -sthe * thedot; - auto stheddot = cthedot * thedot + (cthe * theddot); - auto ctheddot = -(sthedot * thedot) - (sthe * theddot); - auto rotMat = std::make_shared(3, 3); - auto row0 = rotMat->at(0); - row0->atiput(0, ctheddot); - row0->atiput(1, -stheddot); - row0->atiput(2, 0.0); - auto row1 = rotMat->at(1); - row1->atiput(0, stheddot); - row1->atiput(1, ctheddot); - row1->atiput(2, 0.0); - auto row2 = rotMat->at(2); - row2->atiput(0, 0.0); - row2->atiput(1, 0.0); - row2->atiput(2, 0.0); - return rotMat; - } - std::shared_ptr FullMatrixDouble::identitysptr(int n) - { - auto mat = std::make_shared(n, n); - mat->identity(); - return mat; - } - std::shared_ptr FullMatrixFullMatrixDouble::identitysptr(int n) - { - auto mat = std::make_shared(n, n); - mat->identity(); - return mat; - } - std::shared_ptr FullMatrixDouble::tildeMatrix(FColDsptr col) - { - //"tildeMatrix is skew symmetric matrix related to angular velocity and cross product." - if (col->size() != 3) throw std::runtime_error("Column is not of dimension 3"); - auto tilde = std::make_shared(3, 3); - auto c0 = col->at(0); - auto c1 = col->at(1); - auto c2 = col->at(2); - tilde->atijput(0, 0, 0.0); - tilde->atijput(1, 1, 0.0); - tilde->atijput(2, 2, 0.0); - tilde->atijput(1, 2, -c0); - tilde->atijput(0, 2, c1); - tilde->atijput(0, 1, -c2); - tilde->atijput(1, 0, c2); - tilde->atijput(2, 0, -c1); - tilde->atijput(2, 1, c0); - return tilde; - } - void FullMatrixDouble::zeroSelf() - { - for (int i = 0; i < this->size(); i++) { - this->at(i)->zeroSelf(); - } - } - void FullMatrixFullMatrixDouble::zeroSelf() - { - for (int i = 0; i < this->size(); i++) { - this->at(i)->zeroSelf(); - } - } - void FullMatrixFullColumnDouble::zeroSelf() - { - for (int i = 0; i < this->size(); i++) { - this->at(i)->zeroSelf(); - } - } - void FullMatrixDouble::identity() { - this->zeroSelf(); - for (int i = 0; i < this->size(); i++) { - this->at(i)->at(i) = 1.0; - } - } - void FullMatrixFullMatrixDouble::identity() { - assert(false); -// this->zeroSelf(); -// for (int i = 0; i < this->size(); i++) { -// this->at(i)->at(i) = 1.0; -// } - } - // TODO: should there be a FullMatrixFullColumnDouble version of this? - FColsptr FullMatrixDouble::column(int j) { - int n = (int)this->size(); - auto answer = std::make_shared>(n); - for (int i = 0; i < n; i++) { - answer->at(i) = this->at(i)->at(j); - } - return answer; - } - std::shared_ptr FullMatrixDouble::timesFullMatrix(std::shared_ptr fullMat) - { - int m = this->nrow(); - auto answer = std::make_shared(m); - for (int i = 0; i < m; i++) { - answer->at(i) = this->at(i)->timesFullMatrix(fullMat); - } - return answer; - } - std::shared_ptr FullMatrixDouble::timesTransposeFullMatrix(std::shared_ptr fullMat) - { - int nrow = this->nrow(); - auto answer = std::make_shared(nrow); - for (int i = 0; i < nrow; i++) { - answer->at(i) = this->at(i)->timesTransposeFullMatrix(fullMat); - } - return answer; - } -// std::shared_ptr FullMatrixFullMatrixDouble::timesTransposeFullMatrix(std::shared_ptr fullMat) -// { -// int nrow = this->nrow(); -// auto answer = std::make_shared(nrow); -// for (int i = 0; i < nrow; i++) { -// answer->at(i) = this->at(i)->timesTransposeFullMatrixForFMFMDsptr(fullMat); -// } -// return answer; -// } - std::shared_ptr FullMatrixDouble::times(double a) - { - int m = this->nrow(); - auto answer = std::make_shared(m); - for (int i = 0; i < m; i++) { - // auto x = this->at(i); - answer->at(i) = this->at(i)->times(a); - } - return answer; - } - std::shared_ptr FullMatrixFullMatrixDouble::times(double a) - { - // TODO: correct action? - assert(false); - return std::make_shared(); - } - std::shared_ptr FullMatrixFullColumnDouble::times(double a) - { - // TODO: correct action? - assert(false); - return std::make_shared(); - } - std::shared_ptr FullMatrixDouble::transposeTimesFullMatrix(std::shared_ptr fullMat) - { - return this->transpose()->timesFullMatrix(fullMat); - } - std::shared_ptr FullMatrixDouble::plusFullMatrix(std::shared_ptr fullMat) - { - int n = (int)this->size(); - auto answer = std::make_shared(n); - for (int i = 0; i < n; i++) { - answer->at(i) = this->at(i)->plusFullRow(fullMat->at(i)); - } - return answer; - } - std::shared_ptr FullMatrixDouble::minusFullMatrix(std::shared_ptr fullMat) - { - int n = (int)this->size(); - auto answer = std::make_shared(n); - for (int i = 0; i < n; i++) { - answer->at(i) = this->at(i)->minusFullRow(fullMat->at(i)); - } - return answer; - } - std::shared_ptr FullMatrixDouble::transpose() { - int nrow = this->nrow(); - auto ncol = this->ncol(); - auto answer = std::make_shared(ncol, nrow); - for (int i = 0; i < nrow; i++) { - auto& row = this->at(i); - for (int j = 0; j < ncol; j++) { - answer->at(j)->at(i) = row->at(j); - } - } - return answer; - } - std::shared_ptr FullMatrixDouble::negated() - { - return this->times(-1.0); - } - void FullMatrixDouble::symLowerWithUpper() - { - int n = (int)this->size(); - for (int i = 0; i < n; i++) { - for (int j = i + 1; j < n; j++) { - this->at(j)->at(i) = this->at(i)->at(j); - } - } - } - void FullMatrixFullColumnDouble::symLowerWithUpper() - { - int n = (int)this->size(); - for (int i = 0; i < n; i++) { - for (int j = i + 1; j < n; j++) { - this->at(j)->at(i) = this->at(i)->at(j); - } - } - } - void FullMatrixDouble::atiput(int i, FRowsptr fullRow) - { - this->at(i) = fullRow; - } - void FullMatrixDouble::atijput(int i, int j, double value) - { - this->at(i)->atiput(j, value); - } - void FullMatrixDouble::atijputFullColumn(int i1, int j1, FColsptr fullCol) - { - for (int ii = 0; ii < fullCol->size(); ii++) - { - this->at(i1 + ii)->at(j1) = fullCol->at(ii); - } - } - void FullMatrixDouble::atijplusFullRow(int i, int j, FRowsptr fullRow) - { - this->at(i)->atiplusFullRow(j, fullRow); - } - void FullMatrixDouble::atijplusNumber(int i, int j, double value) - { - auto rowi = this->at(i); - rowi->at(j) += value; - } - void FullMatrixDouble::atijminusNumber(int i, int j, double value) - { - auto rowi = this->at(i); - rowi->at(j) -= value; - } - double FullMatrixDouble::sumOfSquares() - { - double sum = 0.0; - for (int i = 0; i < this->size(); i++) - { - sum += this->at(i)->sumOfSquares(); - } - return sum; - } - double FullMatrixFullMatrixDouble::sumOfSquares() - { - double sum = 0.0; - for (int i = 0; i < this->size(); i++) - { - sum += this->at(i)->sumOfSquares(); - } - return sum; - } - double FullMatrixFullColumnDouble::sumOfSquares() - { - double sum = 0.0; - for (int i = 0; i < this->size(); i++) - { - sum += this->at(i)->sumOfSquares(); - } - return sum; - } - std::shared_ptr FullMatrixDouble::copy() - { - auto m = (int)this->size(); - auto answer = std::make_shared(m); - for (int i = 0; i < m; i++) - { - answer->at(i) = this->at(i)->copy(); - } - return answer; - } - FullMatrixDouble FullMatrixDouble::operator+(const FullMatrixDouble fullMat) - { - int n = (int)this->size(); - FullMatrixDouble answer(n); - for (int i = 0; i < n; i++) { - answer.at(i) = this->at(i)->plusFullRow(fullMat.at(i)); - } - return answer; - } - FColsptr FullMatrixDouble::transposeTimesFullColumn(FColsptr fullCol) - { - auto sptr = std::make_shared(*this); - return fullCol->transpose()->timesFullMatrix(sptr)->transpose(); - } - void FullMatrixDouble::magnifySelf(double factor) - { - for (int i = 0; i < this->size(); i++) { - this->at(i)->magnifySelf(factor); - } - } - std::ostream& FullMatrixDouble::printOn(std::ostream& s) const - { - s << "FullMat[" << std::endl; - for (int i = 0; i < this->size(); i++) - { - s << *(this->at(i)) << std::endl; - } - s << "]"; - return s; - } - std::shared_ptr> FullMatrixDouble::asEulerParameters() - { - //"Given [A], compute Euler parameter." +using namespace MbD; - auto traceA = this->trace(); - double dum = 0.0; - double dumSq = 0.0; - //auto qE = CREATE>::With(4); //Cannot use CREATE.h in subclasses of std::vector. Why? - auto qE = std::make_shared>(4); - qE->initialize(); - auto OneMinusTraceDivFour = (1.0 - traceA) / 4.0; - for (int i = 0; i < 3; i++) - { - dumSq = this->at(i)->at(i) / 2.0 + OneMinusTraceDivFour; - dum = (dumSq > 0.0) ? std::sqrt(dumSq) : 0.0; - qE->atiput(i, dum); - } - dumSq = (1.0 + traceA) / 4.0; - dum = (dumSq > 0.0) ? std::sqrt(dumSq) : 0.0; - qE->atiput(3, dum); - double max = 0.0; - int maxE = -1; - for (int i = 0; i < 4; i++) - { - auto num = qE->at(i); - if (max < num) { - max = num; - maxE = i; - } - } - - if (maxE == 0) { - auto FourE = 4.0 * qE->at(0); - qE->atiput(1, (this->at(0)->at(1) + this->at(1)->at(0)) / FourE); - qE->atiput(2, (this->at(0)->at(2) + this->at(2)->at(0)) / FourE); - qE->atiput(3, (this->at(2)->at(1) - this->at(1)->at(2)) / FourE); - } - else if (maxE == 1) { - auto FourE = 4.0 * qE->at(1); - qE->atiput(0, (this->at(0)->at(1) + this->at(1)->at(0)) / FourE); - qE->atiput(2, (this->at(1)->at(2) + this->at(2)->at(1)) / FourE); - qE->atiput(3, (this->at(0)->at(2) - this->at(2)->at(0)) / FourE); - } - else if (maxE == 2) { - auto FourE = 4.0 * qE->at(2); - qE->atiput(0, (this->at(0)->at(2) + this->at(2)->at(0)) / FourE); - qE->atiput(1, (this->at(1)->at(2) + this->at(2)->at(1)) / FourE); - qE->atiput(3, (this->at(1)->at(0) - this->at(0)->at(1)) / FourE); - } - else if (maxE == 3) { - auto FourE = 4.0 * qE->at(3); - qE->atiput(0, (this->at(2)->at(1) - this->at(1)->at(2)) / FourE); - qE->atiput(1, (this->at(0)->at(2) - this->at(2)->at(0)) / FourE); - qE->atiput(2, (this->at(1)->at(0) - this->at(0)->at(1)) / FourE); - } - qE->conditionSelf(); - qE->calc(); - return qE; - } - double FullMatrixDouble::trace() - { - double trace = 0.0; - for (int i = 0; i < this->size(); i++) - { - trace += this->at(i)->at(i); - } - return trace; - } - double FullMatrixDouble::maxMagnitude() - { - double max = 0.0; - for (int i = 0; i < this->size(); i++) - { - double element = this->at(i)->maxMagnitude(); - if (max < element) max = element; - } - return max; - } - double FullMatrixFullMatrixDouble::maxMagnitude() - { - double max = 0.0; - for (int i = 0; i < this->size(); i++) - { - double element = this->at(i)->maxMagnitude(); - if (max < element) max = element; - } - return max; - } - double FullMatrixFullColumnDouble::maxMagnitude() - { - double max = 0.0; - for (int i = 0; i < this->size(); i++) - { - double element = this->at(i)->maxMagnitude(); - if (max < element) max = element; - } - return max; - } - FColsptr FullMatrixDouble::bryantAngles() - { - auto answer = std::make_shared>(3); - auto sthe1y = this->at(0)->at(2); - double the0x, the1y, the2z, cthe0x, sthe0x, y, x; - if (std::abs(sthe1y) > 0.9999) { - if (sthe1y > 0.0) { - the0x = std::atan2(this->at(1)->at(0), this->at(1)->at(1)); - the1y = OS_M_PI / 2.0; - the2z = 0.0; - } - else { - the0x = std::atan2(this->at(2)->at(1), this->at(2)->at(0)); - the1y = OS_M_PI / -2.0; - the2z = 0.0; - } - } - else { - the0x = std::atan2(-this->at(1)->at(2), this->at(2)->at(2)); - cthe0x = std::cos(the0x); - sthe0x = std::sin(the0x); - y = sthe1y; - if (std::abs(cthe0x) > std::abs(sthe0x)) { - x = this->at(2)->at(2) / cthe0x; - } - else { - x = this->at(1)->at(2) / -sthe0x; - } - the1y = std::atan2(y, x); - the2z = std::atan2(-this->at(0)->at(1), this->at(0)->at(0)); - } - answer->atiput(0, the0x); - answer->atiput(1, the1y); - answer->atiput(2, the2z); - return answer; - } - bool FullMatrixDouble::isDiagonal() - { - auto m = this->nrow(); - auto n = this->ncol(); - if (m != n) return false; - for (int i = 0; i < m; i++) - { - auto rowi = this->at(i); - for (int j = 0; j < n; j++) - { - if (i != j && rowi->at(j) != 0) return false; - } - } - return true; - } - bool FullMatrixDouble::isDiagonalToWithin(double ratio) - { - double maxMag = this->maxMagnitude(); - auto tol = ratio * maxMag; - auto nrow = this->nrow(); - if (nrow == this->ncol()) { - for (int i = 0; i < 3; i++) - { - for (int j = i + 1; j < 3; j++) - { - if (std::abs(this->at(i)->at(j)) > tol) return false; - if (std::abs(this->at(j)->at(i)) > tol) return false; - } - } - return true; - } - else { - return false; - } - } - std::shared_ptr FullMatrixDouble::asDiagonalMatrix() - { - int nrow = this->nrow(); - auto diagMat = std::make_shared(nrow); - for (int i = 0; i < nrow; i++) - { - diagMat->atiput(i, this->at(i)->at(i)); - } - return diagMat; - } - void FullMatrixDouble::conditionSelfWithTol(double tol) - { - for (auto row : *this) { - row->conditionSelfWithTol(tol); - } - } -} diff --git a/OndselSolver/FullMatrix.h b/OndselSolver/FullMatrix.h index c21781c..cc3a842 100644 --- a/OndselSolver/FullMatrix.h +++ b/OndselSolver/FullMatrix.h @@ -11,139 +11,731 @@ #include "corecrt_math_defines.h" #include -#include "FullMatrix.ref.h" -#include "FullColumn.ref.h" -#include "FullRow.ref.h" -#include "DiagonalMatrix.ref.h" -#include "EulerParameters.ref.h" #include "RowTypeMatrix.h" -#include "FullRow.h" // now that refs are resolved, go do the full systems -//#include "FullColumn.h" -//#include "DiagonalMatrix.h" -//#include "EulerParameters.h" +#include "FullColumn.h" +#include "FullRow.h" namespace MbD { - // - // FULL MATRIX DOUBLE - // - class FullMatrixDouble : public RowTypeMatrix> { - public: - FullMatrixDouble() = default; - explicit FullMatrixDouble(int m) : RowTypeMatrix>(m) - { - } - FullMatrixDouble(int m, int n) { - for (int i = 0; i < m; i++) { - auto row = std::make_shared>(n); - this->push_back(row); - } - } - FullMatrixDouble(std::initializer_list> listOfRows) { - for (auto& row : listOfRows) - { - this->push_back(row); - } - } - FullMatrixDouble(std::initializer_list> list2D) { - for (auto& rowList : list2D) - { - auto row = std::make_shared>(rowList); - this->push_back(row); - } - } + template + class FullMatrix; + using FMatDsptr = std::shared_ptr>; + template + using FMatsptr = std::shared_ptr>; + template + class FullColumn; + using FColDsptr = std::shared_ptr>; + template + class FullRow; + template + class EulerParameters; + template + class DiagonalMatrix; - std::shared_ptr times(double a); - std::shared_ptr timesTransposeFullMatrix(std::shared_ptr fullMat); - void identity(); - static std::shared_ptr identitysptr(int n); - double sumOfSquares() override; - std::shared_ptr transposeTimesFullMatrix(std::shared_ptr fullMat); - std::shared_ptr timesFullMatrix(std::shared_ptr fullMat); - std::shared_ptr transpose(); - std::shared_ptr negated(); - std::shared_ptr plusFullMatrix(std::shared_ptr fullMat); - static std::shared_ptr rotatex(double angle); - static std::shared_ptr rotatey(double angle); - static std::shared_ptr rotatez(double angle); - static std::shared_ptr rotatexrotDot(double angle, double angledot); - static std::shared_ptr rotateyrotDot(double angle, double angledot); - static std::shared_ptr rotatezrotDot(double angle, double angledot); - static std::shared_ptr rotatexrotDotrotDDot(double angle, double angleDot, double angleDDot); - static std::shared_ptr rotateyrotDotrotDDot(double angle, double angleDot, double angleDDot); - static std::shared_ptr rotatezrotDotrotDDot(double angle, double angleDot, double angleDDot); - static std::shared_ptr tildeMatrix(FColDsptr col); - void zeroSelf() override; - FColsptr column(int j); + using FMatFColDsptr = std::shared_ptr>; + using FMatFMatDsptr = std::shared_ptr>; + using FColFMatDsptr = std::shared_ptr>; - void atiput(int i, FRowsptr fullRow); - void atijput(int i, int j, double value); - std::shared_ptr copy(); - double maxMagnitude() override; - FullMatrixDouble operator+(const FullMatrixDouble fullMat); - std::shared_ptr minusFullMatrix(std::shared_ptr fullMat); - FColsptr transposeTimesFullColumn(FColsptr fullCol); - void symLowerWithUpper(); - void atijputFullColumn(int i1, int j1, FColsptr fullCol); - void atijplusFullRow(int i, int j, FRowsptr fullRow); - void atijplusNumber(int i, int j, double value); - void atijminusNumber(int i, int j, double value); - void magnifySelf(double factor); - std::shared_ptr> asEulerParameters(); - FColsptr bryantAngles(); - double trace(); - bool isDiagonal(); - bool isDiagonalToWithin(double ratio); - std::shared_ptr asDiagonalMatrix(); - void conditionSelfWithTol(double tol); - std::ostream& printOn(std::ostream& s) const override; - FColsptr timesFullColumn(FColsptr fullCol); - FColsptr timesFullColumn(FullColumn* fullCol); - }; - - // - // FULL MATRIX FULL MATRIX DOUBLE - // - class FullMatrixFullMatrixDouble : public RowTypeMatrix> { - public: - FullMatrixFullMatrixDouble() = default; - explicit FullMatrixFullMatrixDouble(int m) : RowTypeMatrix>(m) + template + class FullMatrix : public RowTypeMatrix> + { + public: + FullMatrix() {} + FullMatrix(int m) : RowTypeMatrix>(m) { } - FullMatrixFullMatrixDouble(int m, int n) { - for (int i = 0; i < m; i++) { - auto row = std::make_shared>(n); - this->push_back(row); - } - } - - double maxMagnitude() override; - void zeroSelf() override; - std::shared_ptr times(double a); - // std::shared_ptr timesTransposeFullMatrix(std::shared_ptr fullMat); - double sumOfSquares() override; - void identity(); - static std::shared_ptr identitysptr(int n); - }; - // - // FULL MATRIX FULL COLUMN DOUBLE - // - class FullMatrixFullColumnDouble : public RowTypeMatrix> { - public: - FullMatrixFullColumnDouble() = default; - explicit FullMatrixFullColumnDouble(int m) : RowTypeMatrix>(m) - { + FullMatrix(int m, int n) { + for (int i = 0; i < m; i++) { + auto row = std::make_shared>(n); + this->push_back(row); + } } - FullMatrixFullColumnDouble(int m, int n) { - for (int i = 0; i < m; i++) { - auto row = std::make_shared>(n); - this->push_back(row); - } - } + FullMatrix(std::initializer_list> listOfRows) { + for (auto& row : listOfRows) + { + this->push_back(row); + } + } + FullMatrix(std::initializer_list> list2D) { + for (auto& rowList : list2D) + { + auto row = std::make_shared>(rowList); + this->push_back(row); + } + } + static FMatsptr rotatex(T angle); + static FMatsptr rotatey(T angle); + static FMatsptr rotatez(T angle); + static FMatsptr rotatexrotDot(T angle, T angledot); + static FMatsptr rotateyrotDot(T angle, T angledot); + static FMatsptr rotatezrotDot(T angle, T angledot); + static FMatsptr rotatexrotDotrotDDot(T angle, T angleDot, T angleDDot); + static FMatsptr rotateyrotDotrotDDot(T angle, T angleDot, T angleDDot); + static FMatsptr rotatezrotDotrotDDot(T angle, T angleDot, T angleDDot); + static FMatsptr identitysptr(int n); + static FMatsptr tildeMatrix(FColDsptr col); + void identity(); + FColsptr column(int j); + FColsptr timesFullColumn(FColsptr fullCol); + FColsptr timesFullColumn(FullColumn* fullCol); + FMatsptr timesFullMatrix(FMatsptr fullMat); + FMatsptr timesTransposeFullMatrix(FMatsptr fullMat); + FMatsptr times(T a); + FMatsptr transposeTimesFullMatrix(FMatsptr fullMat); + FMatsptr plusFullMatrix(FMatsptr fullMat); + FMatsptr minusFullMatrix(FMatsptr fullMat); + FMatsptr transpose(); + FMatsptr negated(); + void symLowerWithUpper(); + void atiput(int i, FRowsptr fullRow); + void atijput(int i, int j, T value); + void atijputFullColumn(int i, int j, FColsptr fullCol); + void atijplusFullRow(int i, int j, FRowsptr fullRow); + void atijplusNumber(int i, int j, T value); + void atijminusNumber(int i, int j, T value); + double sumOfSquares() override; + void zeroSelf() override; + FMatsptr copy(); + FullMatrix operator+(const FullMatrix fullMat); + FColsptr transposeTimesFullColumn(const FColsptr fullCol); + void magnifySelf(T factor); + std::shared_ptr> asEulerParameters(); + T trace(); + double maxMagnitude() override; + FColsptr bryantAngles(); + bool isDiagonal(); + bool isDiagonalToWithin(double ratio); + std::shared_ptr> asDiagonalMatrix(); + void conditionSelfWithTol(double tol); - double maxMagnitude() override; - void zeroSelf() override; - double sumOfSquares() override; - void symLowerWithUpper(); - std::shared_ptr times(double a); - }; + std::ostream& printOn(std::ostream& s) const override; + }; + template + inline FMatsptr FullMatrix::rotatex(T the) + { + auto sthe = std::sin(the); + auto cthe = std::cos(the); + auto rotMat = std::make_shared>(3, 3); + auto row0 = rotMat->at(0); + row0->atiput(0, 1.0); + row0->atiput(1, 0.0); + row0->atiput(2, 0.0); + auto row1 = rotMat->at(1); + row1->atiput(0, 0.0); + row1->atiput(1, cthe); + row1->atiput(2, -sthe); + auto row2 = rotMat->at(2); + row2->atiput(0, 0.0); + row2->atiput(1, sthe); + row2->atiput(2, cthe); + return rotMat; + } + template + inline FMatsptr FullMatrix::rotatey(T the) + { + auto sthe = std::sin(the); + auto cthe = std::cos(the); + auto rotMat = std::make_shared>(3, 3); + auto row0 = rotMat->at(0); + row0->atiput(0, cthe); + row0->atiput(1, 0.0); + row0->atiput(2, sthe); + auto row1 = rotMat->at(1); + row1->atiput(0, 0.0); + row1->atiput(1, 1.0); + row1->atiput(2, 0.0); + auto row2 = rotMat->at(2); + row2->atiput(0, -sthe); + row2->atiput(1, 0.0); + row2->atiput(2, cthe); + return rotMat; + } + template + inline FMatsptr FullMatrix::rotatez(T the) + { + auto sthe = std::sin(the); + auto cthe = std::cos(the); + auto rotMat = std::make_shared>(3, 3); + auto row0 = rotMat->at(0); + row0->atiput(0, cthe); + row0->atiput(1, -sthe); + row0->atiput(2, 0.0); + auto row1 = rotMat->at(1); + row1->atiput(0, sthe); + row1->atiput(1, cthe); + row1->atiput(2, 0.0); + auto row2 = rotMat->at(2); + row2->atiput(0, 0.0); + row2->atiput(1, 0.0); + row2->atiput(2, 1.0); + return rotMat; + } + template + inline FMatsptr FullMatrix::rotatexrotDot(T the, T thedot) + { + auto sthe = std::sin(the); + auto cthe = std::cos(the); + auto sthedot = cthe * thedot; + auto cthedot = -sthe * thedot; + auto rotMat = std::make_shared>(3, 3); + auto row0 = rotMat->at(0); + row0->atiput(0, 0.0); + row0->atiput(1, 0.0); + row0->atiput(2, 0.0); + auto row1 = rotMat->at(1); + row1->atiput(0, 0.0); + row1->atiput(1, cthedot); + row1->atiput(2, -sthedot); + auto row2 = rotMat->at(2); + row2->atiput(0, 0.0); + row2->atiput(1, sthedot); + row2->atiput(2, cthedot); + return rotMat; + } + template + inline FMatsptr FullMatrix::rotateyrotDot(T the, T thedot) + { + auto sthe = std::sin(the); + auto cthe = std::cos(the); + auto sthedot = cthe * thedot; + auto cthedot = -sthe * thedot; + auto rotMat = std::make_shared>(3, 3); + auto row0 = rotMat->at(0); + row0->atiput(0, cthedot); + row0->atiput(1, 0.0); + row0->atiput(2, sthedot); + auto row1 = rotMat->at(1); + row1->atiput(0, 0.0); + row1->atiput(1, 0.0); + row1->atiput(2, 0.0); + auto row2 = rotMat->at(2); + row2->atiput(0, -sthedot); + row2->atiput(1, 0.0); + row2->atiput(2, cthedot); + return rotMat; + } + template + inline FMatsptr FullMatrix::rotatezrotDot(T the, T thedot) + { + auto sthe = std::sin(the); + auto cthe = std::cos(the); + auto sthedot = cthe * thedot; + auto cthedot = -sthe * thedot; + auto rotMat = std::make_shared>(3, 3); + auto row0 = rotMat->at(0); + row0->atiput(0, cthedot); + row0->atiput(1, -sthedot); + row0->atiput(2, 0.0); + auto row1 = rotMat->at(1); + row1->atiput(0, sthedot); + row1->atiput(1, cthedot); + row1->atiput(2, 0.0); + auto row2 = rotMat->at(2); + row2->atiput(0, 0.0); + row2->atiput(1, 0.0); + row2->atiput(2, 0.0); + return rotMat; + } + template + inline FMatsptr FullMatrix::rotatexrotDotrotDDot(T the, T thedot, T theddot) + { + auto sthe = std::sin(the); + auto cthe = std::cos(the); + auto sthedot = cthe * thedot; + auto cthedot = -sthe * thedot; + auto stheddot = cthedot * thedot + (cthe * theddot); + auto ctheddot = -(sthedot * thedot) - (sthe * theddot); + auto rotMat = std::make_shared>(3, 3); + auto row0 = rotMat->at(0); + row0->atiput(0, 0.0); + row0->atiput(1, 0.0); + row0->atiput(2, 0.0); + auto row1 = rotMat->at(1); + row1->atiput(0, 0.0); + row1->atiput(1, ctheddot); + row1->atiput(2, -stheddot); + auto row2 = rotMat->at(2); + row2->atiput(0, 0.0); + row2->atiput(1, stheddot); + row2->atiput(2, ctheddot); + return rotMat; + } + template + inline FMatsptr FullMatrix::rotateyrotDotrotDDot(T the, T thedot, T theddot) + { + auto sthe = std::sin(the); + auto cthe = std::cos(the); + auto sthedot = cthe * thedot; + auto cthedot = -sthe * thedot; + auto stheddot = cthedot * thedot + (cthe * theddot); + auto ctheddot = -(sthedot * thedot) - (sthe * theddot); + auto rotMat = std::make_shared>(3, 3); + auto row0 = rotMat->at(0); + row0->atiput(0, ctheddot); + row0->atiput(1, 0.0); + row0->atiput(2, stheddot); + auto row1 = rotMat->at(1); + row1->atiput(0, 0.0); + row1->atiput(1, 0.0); + row1->atiput(2, 0.0); + auto row2 = rotMat->at(2); + row2->atiput(0, -stheddot); + row2->atiput(1, 0.0); + row2->atiput(2, ctheddot); + return rotMat; + } + template + inline FMatsptr FullMatrix::rotatezrotDotrotDDot(T the, T thedot, T theddot) + { + auto sthe = std::sin(the); + auto cthe = std::cos(the); + auto sthedot = cthe * thedot; + auto cthedot = -sthe * thedot; + auto stheddot = cthedot * thedot + (cthe * theddot); + auto ctheddot = -(sthedot * thedot) - (sthe * theddot); + auto rotMat = std::make_shared>(3, 3); + auto row0 = rotMat->at(0); + row0->atiput(0, ctheddot); + row0->atiput(1, -stheddot); + row0->atiput(2, 0.0); + auto row1 = rotMat->at(1); + row1->atiput(0, stheddot); + row1->atiput(1, ctheddot); + row1->atiput(2, 0.0); + auto row2 = rotMat->at(2); + row2->atiput(0, 0.0); + row2->atiput(1, 0.0); + row2->atiput(2, 0.0); + return rotMat; + } + template + inline FMatsptr FullMatrix::identitysptr(int n) + { + auto mat = std::make_shared>(n, n); + mat->identity(); + return mat; + } + template + inline FMatsptr FullMatrix::tildeMatrix(FColDsptr col) + { + //"tildeMatrix is skew symmetric matrix related to angular velocity and cross product." + if (col->size() != 3) throw std::runtime_error("Column is not of dimension 3"); + auto tilde = std::make_shared>(3, 3); + auto c0 = col->at(0); + auto c1 = col->at(1); + auto c2 = col->at(2); + tilde->atijput(0, 0, 0.0); + tilde->atijput(1, 1, 0.0); + tilde->atijput(2, 2, 0.0); + tilde->atijput(1, 2, -c0); + tilde->atijput(0, 2, c1); + tilde->atijput(0, 1, -c2); + tilde->atijput(1, 0, c2); + tilde->atijput(2, 0, -c1); + tilde->atijput(2, 1, c0); + return tilde; + } + template<> + inline void FullMatrix::zeroSelf() + { + for (int i = 0; i < this->size(); i++) { + this->at(i)->zeroSelf(); + } + } + template<> + inline void FullMatrix::identity() { + this->zeroSelf(); + for (int i = 0; i < this->size(); i++) { + this->at(i)->at(i) = 1.0; + } + } + template + inline FColsptr FullMatrix::column(int j) { + int n = (int)this->size(); + auto answer = std::make_shared>(n); + for (int i = 0; i < n; i++) { + answer->at(i) = this->at(i)->at(j); + } + return answer; + } + template + inline FMatsptr FullMatrix::timesFullMatrix(FMatsptr fullMat) + { + int m = this->nrow(); + auto answer = std::make_shared>(m); + for (int i = 0; i < m; i++) { + answer->at(i) = this->at(i)->timesFullMatrix(fullMat); + } + return answer; + } + template + inline FMatsptr FullMatrix::timesTransposeFullMatrix(FMatsptr fullMat) + { + int nrow = this->nrow(); + auto answer = std::make_shared>(nrow); + for (int i = 0; i < nrow; i++) { + answer->at(i) = this->at(i)->timesTransposeFullMatrix(fullMat); + } + return answer; + } + template<> + inline FMatDsptr FullMatrix::times(double a) + { + int m = this->nrow(); + auto answer = std::make_shared>(m); + for (int i = 0; i < m; i++) { + answer->at(i) = this->at(i)->times(a); + } + return answer; + } + template + inline FMatsptr FullMatrix::times(T a) + { + assert(false); + } + template + inline FMatsptr FullMatrix::transposeTimesFullMatrix(FMatsptr fullMat) + { + return this->transpose()->timesFullMatrix(fullMat); + } + template + inline FMatsptr FullMatrix::plusFullMatrix(FMatsptr fullMat) + { + int n = (int)this->size(); + auto answer = std::make_shared>(n); + for (int i = 0; i < n; i++) { + answer->at(i) = this->at(i)->plusFullRow(fullMat->at(i)); + } + return answer; + } + template + inline FMatsptr FullMatrix::minusFullMatrix(FMatsptr fullMat) + { + int n = (int)this->size(); + auto answer = std::make_shared>(n); + for (int i = 0; i < n; i++) { + answer->at(i) = this->at(i)->minusFullRow(fullMat->at(i)); + } + return answer; + } + template + inline FMatsptr FullMatrix::transpose() + { + int nrow = this->nrow(); + auto ncol = this->ncol(); + auto answer = std::make_shared>(ncol, nrow); + for (int i = 0; i < nrow; i++) { + auto& row = this->at(i); + for (int j = 0; j < ncol; j++) { + answer->at(j)->at(i) = row->at(j); + } + } + return answer; + } + template + inline FMatsptr FullMatrix::negated() + { + return this->times(-1.0); + } + template + inline void FullMatrix::symLowerWithUpper() + { + int n = (int)this->size(); + for (int i = 0; i < n; i++) { + for (int j = i + 1; j < n; j++) { + this->at(j)->at(i) = this->at(i)->at(j); + } + } + } + template + inline void FullMatrix::atiput(int i, FRowsptr fullRow) + { + this->at(i) = fullRow; + } + template + inline void FullMatrix::atijput(int i, int j, T value) + { + this->at(i)->atiput(j, value); + } + template + inline void FullMatrix::atijputFullColumn(int i1, int j1, FColsptr fullCol) + { + for (int ii = 0; ii < fullCol->size(); ii++) + { + this->at(i1 + ii)->at(j1) = fullCol->at(ii); + } + } + template + inline void FullMatrix::atijplusFullRow(int i, int j, FRowsptr fullRow) + { + this->at(i)->atiplusFullRow(j, fullRow); + } + template + inline void FullMatrix::atijplusNumber(int i, int j, T value) + { + auto rowi = this->at(i); + rowi->at(j) += value; + } + template + inline void FullMatrix::atijminusNumber(int i, int j, T value) + { + auto rowi = this->at(i); + rowi->at(j) -= value; + } + template<> + inline double FullMatrix::sumOfSquares() + { + double sum = 0.0; + for (int i = 0; i < this->size(); i++) + { + sum += this->at(i)->sumOfSquares(); + } + return sum; + } + template + inline double FullMatrix::sumOfSquares() + { + assert(false); + return 0.0; + } + template + inline void FullMatrix::zeroSelf() + { + assert(false); + } + template + inline FMatsptr FullMatrix::copy() + { + auto m = (int)this->size(); + auto answer = std::make_shared>(m); + for (int i = 0; i < m; i++) + { + answer->at(i) = this->at(i)->copy(); + } + return answer; + } + template + inline FullMatrix FullMatrix::operator+(const FullMatrix fullMat) + { + int n = (int)this->size(); + auto answer = FullMatrix(n); + for (int i = 0; i < n; i++) { + answer.at(i) = this->at(i)->plusFullRow(fullMat.at(i)); + } + return answer; + } + template + inline FColsptr FullMatrix::transposeTimesFullColumn(FColsptr fullCol) + { + auto sptr = std::make_shared>(*this); + return fullCol->transpose()->timesFullMatrix(sptr)->transpose(); + } + template + inline void FullMatrix::magnifySelf(T factor) + { + for (int i = 0; i < this->size(); i++) { + this->at(i)->magnifySelf(factor); + } + } + template + inline std::ostream& FullMatrix::printOn(std::ostream& s) const + { + s << "FullMat[" << std::endl; + for (int i = 0; i < this->size(); i++) + { + s << *(this->at(i)) << std::endl; + } + s << "]"; + return s; + } + template + inline std::shared_ptr> FullMatrix::asEulerParameters() + { + //"Given [A], compute Euler parameter." + + auto traceA = this->trace(); + T dum = 0.0; + T dumSq = 0.0; + //auto qE = CREATE>::With(4); //Cannot use CREATE.h in subclasses of std::vector. Why? + auto qE = std::make_shared>(4); + qE->initialize(); + auto OneMinusTraceDivFour = (1.0 - traceA) / 4.0; + for (int i = 0; i < 3; i++) + { + dumSq = this->at(i)->at(i) / 2.0 + OneMinusTraceDivFour; + dum = (dumSq > 0.0) ? std::sqrt(dumSq) : 0.0; + qE->atiput(i, dum); + } + dumSq = (1.0 + traceA) / 4.0; + dum = (dumSq > 0.0) ? std::sqrt(dumSq) : 0.0; + qE->atiput(3, dum); + T max = 0.0; + int maxE = -1; + for (int i = 0; i < 4; i++) + { + auto num = qE->at(i); + if (max < num) { + max = num; + maxE = i; + } + } + + if (maxE == 0) { + auto FourE = 4.0 * qE->at(0); + qE->atiput(1, (this->at(0)->at(1) + this->at(1)->at(0)) / FourE); + qE->atiput(2, (this->at(0)->at(2) + this->at(2)->at(0)) / FourE); + qE->atiput(3, (this->at(2)->at(1) - this->at(1)->at(2)) / FourE); + } + else if (maxE == 1) { + auto FourE = 4.0 * qE->at(1); + qE->atiput(0, (this->at(0)->at(1) + this->at(1)->at(0)) / FourE); + qE->atiput(2, (this->at(1)->at(2) + this->at(2)->at(1)) / FourE); + qE->atiput(3, (this->at(0)->at(2) - this->at(2)->at(0)) / FourE); + } + else if (maxE == 2) { + auto FourE = 4.0 * qE->at(2); + qE->atiput(0, (this->at(0)->at(2) + this->at(2)->at(0)) / FourE); + qE->atiput(1, (this->at(1)->at(2) + this->at(2)->at(1)) / FourE); + qE->atiput(3, (this->at(1)->at(0) - this->at(0)->at(1)) / FourE); + } + else if (maxE == 3) { + auto FourE = 4.0 * qE->at(3); + qE->atiput(0, (this->at(2)->at(1) - this->at(1)->at(2)) / FourE); + qE->atiput(1, (this->at(0)->at(2) - this->at(2)->at(0)) / FourE); + qE->atiput(2, (this->at(1)->at(0) - this->at(0)->at(1)) / FourE); + } + qE->conditionSelf(); + qE->calc(); + return qE; + } + template + inline T FullMatrix::trace() + { + T trace = 0.0; + for (int i = 0; i < this->size(); i++) + { + trace += this->at(i)->at(i); + } + return trace; + } + template + inline double FullMatrix::maxMagnitude() + { + double max = 0.0; + for (int i = 0; i < this->size(); i++) + { + double element = this->at(i)->maxMagnitude(); + if (max < element) max = element; + } + return max; + } + template + inline FColsptr FullMatrix::bryantAngles() + { + auto answer = std::make_shared>(3); + auto sthe1y = this->at(0)->at(2); + T the0x, the1y, the2z, cthe0x, sthe0x, y, x; + if (std::abs(sthe1y) > 0.9999) { + if (sthe1y > 0.0) { + the0x = std::atan2(this->at(1)->at(0), this->at(1)->at(1)); + the1y = OS_M_PI / 2.0; + the2z = 0.0; + } + else { + the0x = std::atan2(this->at(2)->at(1), this->at(2)->at(0)); + the1y = OS_M_PI / -2.0; + the2z = 0.0; + } + } + else { + the0x = std::atan2(-this->at(1)->at(2), this->at(2)->at(2)); + cthe0x = std::cos(the0x); + sthe0x = std::sin(the0x); + y = sthe1y; + if (std::abs(cthe0x) > std::abs(sthe0x)) { + x = this->at(2)->at(2) / cthe0x; + } + else { + x = this->at(1)->at(2) / -sthe0x; + } + the1y = std::atan2(y, x); + the2z = std::atan2(-this->at(0)->at(1), this->at(0)->at(0)); + } + answer->atiput(0, the0x); + answer->atiput(1, the1y); + answer->atiput(2, the2z); + return answer; + } + template + inline bool FullMatrix::isDiagonal() + { + auto m = this->nrow(); + auto n = this->ncol(); + if (m != n) return false; + for (int i = 0; i < m; i++) + { + auto rowi = this->at(i); + for (int j = 0; j < n; j++) + { + if (i != j && rowi->at(j) != 0) return false; + } + } + return true; + } + template + inline bool FullMatrix::isDiagonalToWithin(double ratio) + { + double maxMag = this->maxMagnitude(); + auto tol = ratio * maxMag; + auto nrow = this->nrow(); + if (nrow == this->ncol()) { + for (int i = 0; i < 3; i++) + { + for (int j = i + 1; j < 3; j++) + { + if (std::abs(this->at(i)->at(j)) > tol) return false; + if (std::abs(this->at(j)->at(i)) > tol) return false; + } + } + return true; + } + else { + return false; + } + } + template + inline std::shared_ptr> FullMatrix::asDiagonalMatrix() + { + int nrow = this->nrow(); + auto diagMat = std::make_shared>(nrow); + for (int i = 0; i < nrow; i++) + { + diagMat->atiput(i, this->at(i)->at(i)); + } + return diagMat; + } + template + inline void FullMatrix::conditionSelfWithTol(double tol) + { + for (auto row : *this) { + row->conditionSelfWithTol(tol); + } + } + template + inline FColsptr FullMatrix::timesFullColumn(FColsptr fullCol) + { + return this->timesFullColumn(fullCol.get()); + } + template + inline FColsptr FullMatrix::timesFullColumn(FullColumn* fullCol) + { + //"a*b = a(i,j)b(j) sum j." + auto nrow = this->nrow(); + auto answer = std::make_shared>(nrow); + for (int i = 0; i < nrow; i++) + { + answer->at(i) = this->at(i)->timesFullColumn(fullCol); + } + return answer; + } } + diff --git a/OndselSolver/FullMatrix.ref.h b/OndselSolver/FullMatrix.ref.h deleted file mode 100644 index a5c0f9e..0000000 --- a/OndselSolver/FullMatrix.ref.h +++ /dev/null @@ -1,19 +0,0 @@ -#pragma once - -#include "FullColumn.ref.h" -#include - -namespace MbD { - class FullMatrixDouble; - class FullMatrixFullMatrixDouble; - class FullMatrixFullColumnDouble; - - using FMatDsptr = std::shared_ptr; - - using FMatFMatDsptr = std::shared_ptr; - - using FColFMatDsptr = std::shared_ptr>; - using FMatFColDsptr = std::shared_ptr; -} - - diff --git a/OndselSolver/FullRow.cpp b/OndselSolver/FullRow.cpp index 109a936..3458277 100644 --- a/OndselSolver/FullRow.cpp +++ b/OndselSolver/FullRow.cpp @@ -7,54 +7,5 @@ ***************************************************************************/ #include "FullRow.h" -#include "FullMatrix.h" -namespace MbD { - template<> - std::shared_ptr FullRow::transposeTimesFullRow(FRowsptr fullRow) - { - //"a*b = a(i)b(j)" - auto nrow = (int)this->size(); - auto answer = std::make_shared(nrow); - for (int i = 0; i < nrow; i++) - { - answer->atiput(i, fullRow->times(this->at(i))); - } - return answer; - } - //template - //inline FMatDsptr FullRow::transposeTimesFullRow(FRowDsptr fullRow) - //{ - // //"a*b = a(i)b(j)" - // auto nrow = (int)this->size(); - // auto answer = std::make_shared(nrow); - // for (int i = 0; i < nrow; i++) - // { - // answer->atiput(i, fullRow->times(this->at(i))); - // } - // return answer; - //} - - template<> - FRowsptr FullRow::timesTransposeFullMatrix(std::shared_ptr fullMat) - { - //"a*bT = a(1,j)b(k,j)" - int ncol = fullMat->nrow(); - auto answer = std::make_shared>(ncol); - for (int k = 0; k < ncol; k++) { - answer->at(k) = this->dot(fullMat->at(k)); - } - return answer; - } - - template<> - FRowsptr FullRow::timesFullMatrix(std::shared_ptr fullMat) - { - FRowsptr answer = fullMat->at(0)->times(this->at(0)); - for (int j = 1; j < this->size(); j++) - { - answer->equalSelfPlusFullRowTimes(fullMat->at(j), this->at(j)); - } - return answer; - } -} +using namespace MbD; diff --git a/OndselSolver/FullRow.h b/OndselSolver/FullRow.h index 4ddb58a..7a989a6 100644 --- a/OndselSolver/FullRow.h +++ b/OndselSolver/FullRow.h @@ -9,10 +9,23 @@ #pragma once #include "FullVector.h" -#include "FullMatrix.ref.h" -#include "FullRow.ref.h" +//#include "FullColumn.h" namespace MbD { + template + class FullRow; + template + using FRowsptr = std::shared_ptr>; + using FRowDsptr = std::shared_ptr>; + template + class FullMatrix; + template + using FMatsptr = std::shared_ptr>; + template + class FullColumn; + template + using FColsptr = std::shared_ptr>; + using ListFRD = std::initializer_list; template class FullRow : public FullVector @@ -30,25 +43,20 @@ namespace MbD { FRowsptr minusFullRow(FRowsptr fullRow); T timesFullColumn(FColsptr fullCol); T timesFullColumn(FullColumn* fullCol); + FRowsptr timesFullMatrix(FMatsptr fullMat); + FRowsptr timesTransposeFullMatrix(FMatsptr fullMat); void equalSelfPlusFullRowTimes(FRowsptr fullRow, double factor); void equalFullRow(FRowsptr fullRow); FColsptr transpose(); FRowsptr copy(); void atiplusFullRow(int j, FRowsptr fullRow); - FMatDsptr transposeTimesFullRow(FRowDsptr fullRow); - std::shared_ptr> clonesptr(); - //double dot(std::shared_ptr> vec); - //double dot(std::shared_ptr> vec); + FMatsptr transposeTimesFullRow(FRowsptr fullRow); double dot(std::shared_ptr> vec); std::shared_ptr> dot(std::shared_ptr>>> vecvec); - std::ostream& printOn(std::ostream& s) const override; - FRowsptr timesTransposeFullMatrix(std::shared_ptr fullMat); - // FRowsptr> timesTransposeFullMatrixForFMFMDsptr(std::shared_ptr fullMat); - FRowsptr timesFullMatrix(std::shared_ptr fullMat); - }; + }; template<> inline FRowDsptr FullRow::times(double a) @@ -73,7 +81,7 @@ namespace MbD { template inline FRowsptr FullRow::plusFullRow(FRowsptr fullRow) { - int n = (int)this->size(); + int n = (int) this->size(); auto answer = std::make_shared>(n); for (int i = 0; i < n; i++) { answer->at(i) = this->at(i) + fullRow->at(i); @@ -83,7 +91,7 @@ namespace MbD { template inline FRowsptr FullRow::minusFullRow(FRowsptr fullRow) { - int n = (int)this->size(); + int n = (int) this->size(); auto answer = std::make_shared>(n); for (int i = 0; i < n; i++) { answer->at(i) = this->at(i) - fullRow->at(i); @@ -106,6 +114,17 @@ namespace MbD { return answer; } template + inline FRowsptr FullRow::timesTransposeFullMatrix(FMatsptr fullMat) + { + //"a*bT = a(1,j)b(k,j)" + int ncol = fullMat->nrow(); + auto answer = std::make_shared>(ncol); + for (int k = 0; k < ncol; k++) { + answer->at(k) = this->dot(fullMat->at(k)); + } + return answer; + } + template inline void FullRow::equalSelfPlusFullRowTimes(FRowsptr fullRow, double factor) { this->equalSelfPlusFullVectortimes(fullRow, factor); @@ -141,9 +160,16 @@ namespace MbD { } } template - inline std::shared_ptr> FullRow::clonesptr() + inline FMatsptr FullRow::transposeTimesFullRow(FRowsptr fullRow) { - return std::make_shared>(*this); + //"a*b = a(i)b(j)" + auto nrow = (int)this->size(); + auto answer = std::make_shared>(nrow); + for (int i = 0; i < nrow; i++) + { + answer->atiput(i, fullRow->times(this->at(i))); + } + return answer; } template inline double FullRow::dot(std::shared_ptr> vec) @@ -182,7 +208,17 @@ namespace MbD { } s << "}"; return s; - }; - + } + template + inline FRowsptr FullRow::timesFullMatrix(FMatsptr fullMat) + { + FRowsptr answer = fullMat->at(0)->times(this->at(0)); + for (int j = 1; j < (int) this->size(); j++) + { + answer->equalSelfPlusFullRowTimes(fullMat->at(j), this->at(j)); + } + return answer; + //return FRowsptr(); + } } diff --git a/OndselSolver/FullRow.ref.h b/OndselSolver/FullRow.ref.h deleted file mode 100644 index a553b4e..0000000 --- a/OndselSolver/FullRow.ref.h +++ /dev/null @@ -1,15 +0,0 @@ -#pragma once - -#include "FullColumn.ref.h" -#include - -namespace MbD { - template - class FullRow; - - template - using FRowsptr = std::shared_ptr>; - using FRowDsptr = std::shared_ptr>; - - using ListFRD = std::initializer_list; -} \ No newline at end of file diff --git a/OndselSolver/GeneralSpline.cpp b/OndselSolver/GeneralSpline.cpp index f8d717d..0e311ba 100644 --- a/OndselSolver/GeneralSpline.cpp +++ b/OndselSolver/GeneralSpline.cpp @@ -119,7 +119,7 @@ void MbD::GeneralSpline::computeDerivatives() } auto solver = CREATE::With(); auto derivsVector = solver->solvewithsaveOriginal(matrix, bvector, false); - derivs = std::make_shared(n, p); + derivs = std::make_shared>(n, p); auto hmaxpowers = std::make_shared>(p); for (int j = 0; j < p; j++) { diff --git a/OndselSolver/Item.h b/OndselSolver/Item.h index 7b02d56..7fff0c0 100644 --- a/OndselSolver/Item.h +++ b/OndselSolver/Item.h @@ -10,10 +10,6 @@ #include #include -#include "FullColumn.ref.h" -#include "FullRow.ref.h" -#include "DiagonalMatrix.ref.h" -#include "FullMatrix.ref.h" #include "FullColumn.h" #include "FullRow.h" #include "FullMatrix.h" diff --git a/OndselSolver/LDUFullMat.cpp b/OndselSolver/LDUFullMat.cpp index 94f1cb5..5a7a4c7 100644 --- a/OndselSolver/LDUFullMat.cpp +++ b/OndselSolver/LDUFullMat.cpp @@ -116,7 +116,7 @@ FMatDsptr LDUFullMat::inversesaveOriginal(FMatDsptr fullMat, bool saveOriginal) this->decomposesaveOriginal(fullMat, saveOriginal); rightHandSideB = std::make_shared>(m); - auto matrixAinverse = std::make_shared (m, n); + auto matrixAinverse = std::make_shared >(m, n); for (int j = 0; j < n; j++) { rightHandSideB->zeroSelf(); diff --git a/OndselSolver/MBDynBody.cpp b/OndselSolver/MBDynBody.cpp index 9a218a3..40ab5d0 100644 --- a/OndselSolver/MBDynBody.cpp +++ b/OndselSolver/MBDynBody.cpp @@ -60,7 +60,7 @@ void MbD::MBDynBody::readInertiaMatrix(std::vector& args) { auto parser = std::make_shared(); parser->variables = mbdynVariables(); - aJmat = std::make_shared(3, 3); + aJmat = std::make_shared>(3, 3); auto str = args.at(0); //Must copy string if (str.find("diag") != std::string::npos) { args.erase(args.begin()); @@ -89,7 +89,7 @@ void MbD::MBDynBody::createASMT() if (aJmat->isDiagonalToWithin(1.0e-6)) { asmtMassMarker->setMomentOfInertias(aJmat->asDiagonalMatrix()); asmtMassMarker->setPosition3D(rPcmP); - asmtMassMarker->setRotationMatrix(FullMatrixDouble::identitysptr(3)); + asmtMassMarker->setRotationMatrix(FullMatrix::identitysptr(3)); auto asmtPart = asmtAssembly()->partPartialNamed(nodeName); asmtPart->setPrincipalMassMarker(asmtMassMarker); } @@ -98,7 +98,7 @@ void MbD::MBDynBody::createASMT() solver->setm(mass); solver->setJPP(aJmat); solver->setrPoP(rPcmP); - solver->setAPo(FullMatrixDouble::identitysptr(3)); + solver->setAPo(FullMatrix::identitysptr(3)); solver->setrPcmP(rPcmP); solver->calc(); asmtMassMarker->setMomentOfInertias(solver->aJpp); diff --git a/OndselSolver/MBDynItem.cpp b/OndselSolver/MBDynItem.cpp index b73f129..ffa1e0d 100644 --- a/OndselSolver/MBDynItem.cpp +++ b/OndselSolver/MBDynItem.cpp @@ -242,7 +242,7 @@ FColDsptr MbD::MBDynItem::readBasicPosition(std::vector& args) FMatDsptr MbD::MBDynItem::readOrientation(std::vector& args) { - auto aAOf = FullMatrixDouble::identitysptr(3); + auto aAOf = FullMatrix::identitysptr(3); if (args.empty()) return aAOf; auto str = args.at(0); //Must copy string if (str.find("reference") != std::string::npos) { @@ -292,7 +292,7 @@ FMatDsptr MbD::MBDynItem::readBasicOrientation(std::vector& args) } if (str.find("eye") != std::string::npos) { args.erase(args.begin()); - auto aAFf = FullMatrixDouble::identitysptr(3); + auto aAFf = FullMatrix::identitysptr(3); return aAFf; } auto iss = std::istringstream(str); @@ -341,7 +341,7 @@ FMatDsptr MbD::MBDynItem::readBasicOrientation(std::vector& args) else { assert(false); } - auto aAFf = FullMatrixDouble::identitysptr(3); + auto aAFf = FullMatrix::identitysptr(3); aAFf->atijputFullColumn(0, 0, vecX); aAFf->atijputFullColumn(0, 1, vecY); aAFf->atijputFullColumn(0, 2, vecZ); @@ -390,13 +390,13 @@ FMatDsptr MbD::MBDynItem::readBasicOrientation(std::vector& args) else { assert(false); } - auto aAFf = FullMatrixDouble::identitysptr(3); + auto aAFf = FullMatrix::identitysptr(3); aAFf->atijputFullColumn(0, 0, vecX); aAFf->atijputFullColumn(0, 1, vecY); aAFf->atijputFullColumn(0, 2, vecZ); return aAFf; } - auto aAFf = FullMatrixDouble::identitysptr(3); + auto aAFf = FullMatrix::identitysptr(3); for (int i = 0; i < 3; i++) { auto& rowi = aAFf->at(i); diff --git a/OndselSolver/MBDynJoint.cpp b/OndselSolver/MBDynJoint.cpp index 158c696..17a677e 100644 --- a/OndselSolver/MBDynJoint.cpp +++ b/OndselSolver/MBDynJoint.cpp @@ -119,7 +119,7 @@ void MbD::MBDynJoint::parseMBDyn(std::string line) mkr1->owner = this; mkr1->nodeStr = "Assembly"; mkr1->rPmP = std::make_shared>(3); - mkr1->aAPm = FullMatrixDouble::identitysptr(3); + mkr1->aAPm = FullMatrix::identitysptr(3); readClampMarkerJ(arguments); return; } diff --git a/OndselSolver/MBDynMarker.cpp b/OndselSolver/MBDynMarker.cpp index b540a0e..758933e 100644 --- a/OndselSolver/MBDynMarker.cpp +++ b/OndselSolver/MBDynMarker.cpp @@ -11,7 +11,7 @@ using namespace MbD; void MbD::MBDynMarker::parseMBDyn(std::vector& args) { rPmP = std::make_shared>(3); - aAPm = FullMatrixDouble::identitysptr(3); + aAPm = FullMatrix::identitysptr(3); if (args.empty()) return; auto str = args.at(0); //Must copy string if (str.find("reference") != std::string::npos) { @@ -37,7 +37,7 @@ void MbD::MBDynMarker::parseMBDynClamp(std::vector& args) //rOmO = rOPO + aAOP*rPmP //aAOm = aAOP * aAPm auto rOmO = std::make_shared>(3); - auto aAOm = FullMatrixDouble::identitysptr(3); + auto aAOm = FullMatrix::identitysptr(3); auto rOPO = readPosition(args); auto aAOP = readOrientation(args); rPmP = aAOP->transposeTimesFullColumn(rOmO->minusFullColumn(rOPO)); diff --git a/OndselSolver/MBDynStructural.cpp b/OndselSolver/MBDynStructural.cpp index e4ae6fe..ab666d3 100644 --- a/OndselSolver/MBDynStructural.cpp +++ b/OndselSolver/MBDynStructural.cpp @@ -12,7 +12,7 @@ using namespace MbD; MbD::MBDynStructural::MBDynStructural() { rOfO = std::make_shared>(3); - aAOf = FullMatrixDouble::identitysptr(3); + aAOf = FullMatrix::identitysptr(3); vOfO = std::make_shared>(3); omeOfO = std::make_shared>(3); } diff --git a/OndselSolver/MarkerFrame.cpp b/OndselSolver/MarkerFrame.cpp index fb9eafa..febbf6f 100644 --- a/OndselSolver/MarkerFrame.cpp +++ b/OndselSolver/MarkerFrame.cpp @@ -31,7 +31,7 @@ System* MarkerFrame::root() void MarkerFrame::initialize() { - prOmOpE = std::make_shared(3, 4); + prOmOpE = std::make_shared>(3, 4); pAOmpE = std::make_shared>(4); endFrames = std::make_shared>(); auto endFrm = CREATE::With(); diff --git a/OndselSolver/MarkerFrame.h b/OndselSolver/MarkerFrame.h index 6950310..5ee8cea 100644 --- a/OndselSolver/MarkerFrame.h +++ b/OndselSolver/MarkerFrame.h @@ -72,9 +72,9 @@ namespace MbD { PartFrame* partFrame; //Use raw pointer when pointing backwards. FColDsptr rpmp = std::make_shared>(3); - FMatDsptr aApm = FullMatrixDouble::identitysptr(3); + FMatDsptr aApm = FullMatrix::identitysptr(3); FColDsptr rOmO = std::make_shared>(3); - FMatDsptr aAOm = FullMatrixDouble::identitysptr(3); + FMatDsptr aAOm = FullMatrix::identitysptr(3); FMatDsptr prOmOpE; FColFMatDsptr pAOmpE; FMatFColDsptr pprOmOpEpE; diff --git a/OndselSolver/MomentOfInertiaSolver.cpp b/OndselSolver/MomentOfInertiaSolver.cpp index 883d27d..91c2dd9 100644 --- a/OndselSolver/MomentOfInertiaSolver.cpp +++ b/OndselSolver/MomentOfInertiaSolver.cpp @@ -7,7 +7,7 @@ using namespace MbD; void MbD::MomentOfInertiaSolver::example1() { - auto aJpp = std::make_shared(ListListD{ + auto aJpp = std::make_shared>(ListListD{ { 1, 0, 0 }, { 0, 2, 0 }, { 0, 0, 3 } @@ -56,8 +56,8 @@ void MbD::MomentOfInertiaSolver::doFullPivoting(int p) auto mag = std::abs(aij); if (mag > max) { max = mag; - pivotRow = (int)i; - pivotCol = (int)j; + pivotRow = i; + pivotCol = j; } } } @@ -179,10 +179,10 @@ void MbD::MomentOfInertiaSolver::calcJoo() if (!rPoP) { rPoP = rPcmP; - aAPo = FullMatrixDouble::identitysptr(3); + aAPo = FullMatrix::identitysptr(3); } - auto rocmPtilde = FullMatrixDouble::tildeMatrix(rPcmP->minusFullColumn(rPoP)); - auto rPoPtilde = FullMatrixDouble::tildeMatrix(rPoP); + auto rocmPtilde = FullMatrix::tildeMatrix(rPcmP->minusFullColumn(rPoP)); + auto rPoPtilde = FullMatrix::tildeMatrix(rPoP); auto term1 = aJPP; auto term21 = rPoPtilde->timesFullMatrix(rPoPtilde); auto term22 = rPoPtilde->timesFullMatrix(rocmPtilde); @@ -197,7 +197,7 @@ void MbD::MomentOfInertiaSolver::calcJpp() { //"aJcmP = aJPP + mass*(rPcmPTilde*rPcmPTilde)" - auto rPcmPtilde = FullMatrixDouble::tildeMatrix(rPcmP); + auto rPcmPtilde = FullMatrix::tildeMatrix(rPcmP); aJcmP = aJPP->plusFullMatrix(rPcmPtilde->timesFullMatrix(rPcmPtilde)->times(m)); aJcmP->symLowerWithUpper(); aJcmP->conditionSelfWithTol(aJcmP->maxMagnitude() * 1.0e-6); @@ -217,7 +217,7 @@ void MbD::MomentOfInertiaSolver::calcAPp() auto lam2 = aJpp->at(2); if (lam0 == lam1) { if (lam1 == lam2) { - aAPp = FullMatrixDouble::identitysptr(3); + aAPp = FullMatrix::identitysptr(3); } else { eigenvector1 = eigenvectorFor(lam1); @@ -238,7 +238,7 @@ void MbD::MomentOfInertiaSolver::calcAPp() if (eigenvector1->at(1) < 0.0) eigenvector1->negateSelf(); eigenvector2 = eigenvector0->cross(eigenvector1); } - aAPp = std::make_shared(3, 3); + aAPp = std::make_shared>(3, 3); aAPp->atijputFullColumn(0, 0, eigenvector0); aAPp->atijputFullColumn(0, 1, eigenvector1); aAPp->atijputFullColumn(0, 2, eigenvector2); @@ -295,7 +295,7 @@ void MbD::MomentOfInertiaSolver::calcJppFromDiagJcmP() //"Eigenvalues are orders from smallest to largest." double average; - auto sortedJ = std::make_shared(); + auto sortedJ = std::make_shared>(); sortedJ->push_back(aJcmP->at(0)->at(0)); sortedJ->push_back(aJcmP->at(1)->at(1)); sortedJ->push_back(aJcmP->at(2)->at(2)); @@ -326,7 +326,7 @@ void MbD::MomentOfInertiaSolver::calcJppFromDiagJcmP() lam2 = average; } } - aJpp = std::make_shared(ListD{ lam0, lam1, lam2 }); + aJpp = std::make_shared>(ListD{ lam0, lam1, lam2 }); } void MbD::MomentOfInertiaSolver::calcJppFromFullJcmP() @@ -351,7 +351,7 @@ void MbD::MomentOfInertiaSolver::calcJppFromFullJcmP() auto phiDiv3 = modifiedArcCos(-q / std::sqrt(-p * p * p)) / 3.0; auto twoSqrtMinusp = 2.0 * std::sqrt(-p); auto piDiv3 = OS_M_PI / 3.0; - auto sortedJ = std::make_shared(); + auto sortedJ = std::make_shared>(); sortedJ->push_back(twoSqrtMinusp * std::cos(phiDiv3)); sortedJ->push_back(twoSqrtMinusp * -std::cos(phiDiv3 + piDiv3)); sortedJ->push_back(twoSqrtMinusp * -std::cos(phiDiv3 - piDiv3)); @@ -382,7 +382,7 @@ void MbD::MomentOfInertiaSolver::calcJppFromFullJcmP() lam2 = average; } } - aJpp = std::make_shared(ListD{ lam0, lam1, lam2 }); + aJpp = std::make_shared>(ListD{ lam0, lam1, lam2 }); } double MbD::MomentOfInertiaSolver::modifiedArcCos(double val) diff --git a/OndselSolver/OndselSolver.vcxproj b/OndselSolver/OndselSolver.vcxproj index a08a59a..89ac539 100644 --- a/OndselSolver/OndselSolver.vcxproj +++ b/OndselSolver/OndselSolver.vcxproj @@ -533,7 +533,6 @@ - @@ -543,9 +542,6 @@ - - - diff --git a/OndselSolver/OndselSolver.vcxproj.filters b/OndselSolver/OndselSolver.vcxproj.filters index 2327bed..271f825 100644 --- a/OndselSolver/OndselSolver.vcxproj.filters +++ b/OndselSolver/OndselSolver.vcxproj.filters @@ -1814,21 +1814,9 @@ Header Files - - Header Files - - - Header Files - - - Header Files - Header Files - - Header Files - Header Files diff --git a/OndselSolver/OrbitAngleZIeqcJec.cpp b/OndselSolver/OrbitAngleZIeqcJec.cpp index a5c3334..1dd5298 100644 --- a/OndselSolver/OrbitAngleZIeqcJec.cpp +++ b/OndselSolver/OrbitAngleZIeqcJec.cpp @@ -138,9 +138,9 @@ void MbD::OrbitAngleZIeqcJec::initialize() OrbitAngleZIecJec::initialize(); pthezpXI = std::make_shared>(3); pthezpEI = std::make_shared>(4); - ppthezpXIpXI = std::make_shared(3, 3); - ppthezpXIpEI = std::make_shared(3, 4); - ppthezpEIpEI = std::make_shared(4, 4); + ppthezpXIpXI = std::make_shared>(3, 3); + ppthezpXIpEI = std::make_shared>(3, 4); + ppthezpEIpEI = std::make_shared>(4, 4); } FMatDsptr MbD::OrbitAngleZIeqcJec::ppvaluepEIpEI() diff --git a/OndselSolver/OrbitAngleZIeqcJeqc.cpp b/OndselSolver/OrbitAngleZIeqcJeqc.cpp index b29ea9e..79709c6 100644 --- a/OndselSolver/OrbitAngleZIeqcJeqc.cpp +++ b/OndselSolver/OrbitAngleZIeqcJeqc.cpp @@ -242,13 +242,13 @@ void MbD::OrbitAngleZIeqcJeqc::initialize() OrbitAngleZIeqcJec::initialize(); pthezpXJ = std::make_shared>(3); pthezpEJ = std::make_shared>(4); - ppthezpXIpXJ = std::make_shared(3, 3); - ppthezpXIpEJ = std::make_shared(3, 4); - ppthezpEIpXJ = std::make_shared(4, 3); - ppthezpEIpEJ = std::make_shared(4, 4); - ppthezpXJpXJ = std::make_shared(3, 3); - ppthezpXJpEJ = std::make_shared(3, 4); - ppthezpEJpEJ = std::make_shared(4, 4); + ppthezpXIpXJ = std::make_shared>(3, 3); + ppthezpXIpEJ = std::make_shared>(3, 4); + ppthezpEIpXJ = std::make_shared>(4, 3); + ppthezpEIpEJ = std::make_shared>(4, 4); + ppthezpXJpXJ = std::make_shared>(3, 3); + ppthezpXJpEJ = std::make_shared>(3, 4); + ppthezpEJpEJ = std::make_shared>(4, 4); } FMatDsptr MbD::OrbitAngleZIeqcJeqc::ppvaluepEIpEJ() diff --git a/OndselSolver/Part.cpp b/OndselSolver/Part.cpp index 0bdcec9..35063e0 100644 --- a/OndselSolver/Part.cpp +++ b/OndselSolver/Part.cpp @@ -35,18 +35,18 @@ void Part::initialize() partFrame = CREATE::With(); partFrame->setPart(this); pTpE = std::make_shared>(4); - ppTpEpE = std::make_shared(4, 4); - ppTpEpEdot = std::make_shared(4, 4); + ppTpEpE = std::make_shared>(4, 4); + ppTpEpEdot = std::make_shared>(4, 4); } void Part::initializeLocally() { partFrame->initializeLocally(); if (m > 0) { - mX = std::make_shared(3, m); + mX = std::make_shared>(3, m); } else { - mX = std::make_shared(3, 0.0); + mX = std::make_shared>(3, 0.0); } } @@ -271,8 +271,8 @@ void Part::fillqsuWeights(DiagMatDsptr diagMat) auto aJiMax = this->root()->maximumMomentOfInertia(); double minw = 1.0e3; double maxw = 1.0e6; - auto wqX = std::make_shared(3); - auto wqE = std::make_shared(4); + auto wqX = std::make_shared>(3); + auto wqE = std::make_shared>(4); if (mMax == 0) { mMax = 1.0; } for (int i = 0; i < 3; i++) { @@ -319,8 +319,8 @@ void Part::fillqsudotWeights(DiagMatDsptr diagMat) if (maxInertia == 0) maxInertia = 1.0; double minw = 1.0e-12 * maxInertia; double maxw = maxInertia; - auto wqXdot = std::make_shared(3); - auto wqEdot = std::make_shared(4); + auto wqXdot = std::make_shared>(3); + auto wqEdot = std::make_shared>(4); for (int i = 0; i < 3; i++) { wqXdot->at(i) = (maxw * m / maxInertia) + minw; diff --git a/OndselSolver/Part.h b/OndselSolver/Part.h index ee6c255..a5a0442 100644 --- a/OndselSolver/Part.h +++ b/OndselSolver/Part.h @@ -11,9 +11,6 @@ #include #include "Item.h" -#include "FullColumn.ref.h" -#include "FullMatrix.ref.h" -#include "DiagonalMatrix.ref.h" #include "EulerParametersDot.h" namespace MbD { diff --git a/OndselSolver/RackPinConstraintIqcJc.cpp b/OndselSolver/RackPinConstraintIqcJc.cpp index 6ca70f3..1ce7c5c 100644 --- a/OndselSolver/RackPinConstraintIqcJc.cpp +++ b/OndselSolver/RackPinConstraintIqcJc.cpp @@ -18,8 +18,8 @@ MbD::RackPinConstraintIqcJc::RackPinConstraintIqcJc(EndFrmsptr frmi, EndFrmsptr { pGpXI = std::make_shared>(3); pGpEI = std::make_shared>(4); - ppGpXIpEI = std::make_shared(3, 4); - ppGpEIpEI = std::make_shared(4, 4); + ppGpXIpEI = std::make_shared>(3, 4); + ppGpEIpEI = std::make_shared>(4, 4); } void MbD::RackPinConstraintIqcJc::initxIeJeIe() diff --git a/OndselSolver/RackPinConstraintIqcJqc.cpp b/OndselSolver/RackPinConstraintIqcJqc.cpp index 76317f4..8291ef5 100644 --- a/OndselSolver/RackPinConstraintIqcJqc.cpp +++ b/OndselSolver/RackPinConstraintIqcJqc.cpp @@ -18,9 +18,9 @@ MbD::RackPinConstraintIqcJqc::RackPinConstraintIqcJqc(EndFrmsptr frmi, EndFrmspt { pGpXJ = std::make_shared>(3); pGpEJ = std::make_shared>(4); - ppGpEIpXJ = std::make_shared(4, 3); - ppGpEIpEJ = std::make_shared(4, 4); - ppGpEJpEJ = std::make_shared(4, 4); + ppGpEIpXJ = std::make_shared>(4, 3); + ppGpEIpEJ = std::make_shared>(4, 4); + ppGpEJpEJ = std::make_shared>(4, 4); } void MbD::RackPinConstraintIqcJqc::initxIeJeIe() diff --git a/OndselSolver/RowTypeMatrix.h b/OndselSolver/RowTypeMatrix.h index f4e2504..39ffa09 100644 --- a/OndselSolver/RowTypeMatrix.h +++ b/OndselSolver/RowTypeMatrix.h @@ -9,7 +9,6 @@ #pragma once #include "Array.h" -#include "FullRow.ref.h" namespace MbD { diff --git a/OndselSolver/ScrewConstraintIqcJc.cpp b/OndselSolver/ScrewConstraintIqcJc.cpp index fde053f..5c39af4 100644 --- a/OndselSolver/ScrewConstraintIqcJc.cpp +++ b/OndselSolver/ScrewConstraintIqcJc.cpp @@ -20,8 +20,8 @@ MbD::ScrewConstraintIqcJc::ScrewConstraintIqcJc(EndFrmsptr frmi, EndFrmsptr frmj { pGpXI = std::make_shared>(3); pGpEI = std::make_shared>(4); - ppGpXIpEI = std::make_shared(3, 4); - ppGpEIpEI = std::make_shared(4, 4); + ppGpXIpEI = std::make_shared>(3, 4); + ppGpEIpEI = std::make_shared>(4, 4); } void MbD::ScrewConstraintIqcJc::initzIeJeIe() diff --git a/OndselSolver/ScrewConstraintIqcJqc.cpp b/OndselSolver/ScrewConstraintIqcJqc.cpp index fc25d6c..bdc6ab7 100644 --- a/OndselSolver/ScrewConstraintIqcJqc.cpp +++ b/OndselSolver/ScrewConstraintIqcJqc.cpp @@ -20,9 +20,9 @@ MbD::ScrewConstraintIqcJqc::ScrewConstraintIqcJqc(EndFrmsptr frmi, EndFrmsptr fr { pGpXJ = std::make_shared>(3); pGpEJ = std::make_shared>(4); - ppGpEIpXJ = std::make_shared(4, 3); - ppGpEIpEJ = std::make_shared(4, 4); - ppGpEJpEJ = std::make_shared(4, 4); + ppGpEIpXJ = std::make_shared>(4, 3); + ppGpEIpEJ = std::make_shared>(4, 4); + ppGpEJpEJ = std::make_shared>(4, 4); } void MbD::ScrewConstraintIqcJqc::initzIeJeIe() diff --git a/OndselSolver/StableBackwardDifference.cpp b/OndselSolver/StableBackwardDifference.cpp index 29c4f43..48a27dd 100644 --- a/OndselSolver/StableBackwardDifference.cpp +++ b/OndselSolver/StableBackwardDifference.cpp @@ -48,7 +48,7 @@ FColDsptr MbD::StableBackwardDifference::derivativepresentpastpresentDerivativep void StableBackwardDifference::instantiateTaylorMatrix() { if (taylorMatrix == nullptr || (taylorMatrix->nrow() != (order))) { - taylorMatrix = std::make_shared(order, order); + taylorMatrix = std::make_shared>(order, order); } } @@ -79,7 +79,7 @@ FColDsptr MbD::StableBackwardDifference::derivativepresentpast(int deriv, FColDs //"Answer ith derivative given present value and past values." if (deriv == 0) { - return y->cloneFcSptr(); + return std::static_pointer_cast>(y->clonesptr()); } else { if (deriv <= order) { diff --git a/OndselSolver/VelICSolver.cpp b/OndselSolver/VelICSolver.cpp index 07f582d..b36ab99 100644 --- a/OndselSolver/VelICSolver.cpp +++ b/OndselSolver/VelICSolver.cpp @@ -59,7 +59,7 @@ void VelICSolver::runBasic() this->assignEquationNumbers(); system->partsJointsMotionsDo([](std::shared_ptr item) { item->useEquationNumbers(); }); auto qsudotOld = std::make_shared>(nqsu); - auto qsudotWeights = std::make_shared(nqsu); + auto qsudotWeights = std::make_shared>(nqsu); errorVector = std::make_shared>(n); jacobian = std::make_shared>(n, n); system->partsJointsMotionsDo([&](std::shared_ptr item) { item->fillqsudot(qsudotOld); }); From 2b7fba1b0efd5931e28b7bf5ff4604d913ce075b Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Thu, 30 Nov 2023 21:48:26 -0700 Subject: [PATCH 02/12] Edits for MacOS --- OndselSolver/FullColumn.cpp | 2 -- OndselSolver/MBDynBody.h | 2 +- OndselSolver/MBDynDrive.h | 2 +- OndselSolver/MBDynGravity.h | 2 +- OndselSolver/MBDynItem.cpp | 7 ++++++- OndselSolver/MBDynItem.h | 1 + OndselSolver/MBDynJoint.h | 2 +- OndselSolver/MBDynReference.h | 2 +- OndselSolver/MBDynStructural.h | 2 +- 9 files changed, 13 insertions(+), 9 deletions(-) diff --git a/OndselSolver/FullColumn.cpp b/OndselSolver/FullColumn.cpp index ae0ee14..e689f73 100644 --- a/OndselSolver/FullColumn.cpp +++ b/OndselSolver/FullColumn.cpp @@ -9,7 +9,5 @@ #include #include "FullColumn.h" -#include "FullRow.h" -#include "FullMatrix.h" using namespace MbD; diff --git a/OndselSolver/MBDynBody.h b/OndselSolver/MBDynBody.h index f93df5a..2151fdc 100644 --- a/OndselSolver/MBDynBody.h +++ b/OndselSolver/MBDynBody.h @@ -14,7 +14,7 @@ namespace MbD { { public: void initialize() override; - void parseMBDyn(std::string line); + void parseMBDyn(std::string line) override; void readMass(std::vector& args); void readInertiaMatrix(std::vector& args); void createASMT() override; diff --git a/OndselSolver/MBDynDrive.h b/OndselSolver/MBDynDrive.h index a4dad0a..31425e4 100644 --- a/OndselSolver/MBDynDrive.h +++ b/OndselSolver/MBDynDrive.h @@ -13,7 +13,7 @@ namespace MbD { class MBDynDrive : public MBDynElement { public: - void parseMBDyn(std::string line); + void parseMBDyn(std::string line) override; void readFunction(std::vector& args); void createASMT() override; diff --git a/OndselSolver/MBDynGravity.h b/OndselSolver/MBDynGravity.h index 4dd9f0e..0ae1a80 100644 --- a/OndselSolver/MBDynGravity.h +++ b/OndselSolver/MBDynGravity.h @@ -13,7 +13,7 @@ namespace MbD { class MBDynGravity : public MBDynElement { public: - void parseMBDyn(std::string line); + void parseMBDyn(std::string line) override; void readFunction(std::vector& args); void createASMT() override; diff --git a/OndselSolver/MBDynItem.cpp b/OndselSolver/MBDynItem.cpp index ffa1e0d..93d822c 100644 --- a/OndselSolver/MBDynItem.cpp +++ b/OndselSolver/MBDynItem.cpp @@ -34,6 +34,11 @@ void MbD::MBDynItem::parseMBDyn(std::vector& lines) assert(false); } +void MbD::MBDynItem::parseMBDyn(std::string line) +{ + assert(false); +} + std::vector MbD::MBDynItem::collectArgumentsFor(std::string title, std::string& statement) { size_t previousPos = 0; @@ -67,7 +72,7 @@ std::vector MbD::MBDynItem::collectArgumentsFor(std::string title, //Need to find matching '"' auto it = std::find_if(arguments.begin() + 1, arguments.end(), [](const std::string& s) { auto nn = std::count(s.begin(), s.end(), '"'); - if ((nn % 2) == 1) return true; + return (nn % 2) == 1; }); std::vector needToCombineArgs(arguments.begin(), it + 1); arguments.erase(arguments.begin(), it + 1); diff --git a/OndselSolver/MBDynItem.h b/OndselSolver/MBDynItem.h index a994c01..c8082ec 100644 --- a/OndselSolver/MBDynItem.h +++ b/OndselSolver/MBDynItem.h @@ -31,6 +31,7 @@ namespace MbD { void noop(); //void setName(std::string str); virtual void parseMBDyn(std::vector& lines); + virtual void parseMBDyn(std::string line); static std::vector collectArgumentsFor(std::string title, std::string& statement); std::vector::iterator findLineWith(std::vector& lines, std::vector& tokens); bool lineHasTokens(const std::string& line, std::vector& tokens); diff --git a/OndselSolver/MBDynJoint.h b/OndselSolver/MBDynJoint.h index be2f68e..ccca9a1 100644 --- a/OndselSolver/MBDynJoint.h +++ b/OndselSolver/MBDynJoint.h @@ -15,7 +15,7 @@ namespace MbD { { public: void initialize() override; - void parseMBDyn(std::string line); + void parseMBDyn(std::string line) override; void readMarkerI(std::vector& args); void readMarkerJ(std::vector& args); void readClampMarkerJ(std::vector& args); diff --git a/OndselSolver/MBDynReference.h b/OndselSolver/MBDynReference.h index 5867816..dfe9f90 100644 --- a/OndselSolver/MBDynReference.h +++ b/OndselSolver/MBDynReference.h @@ -15,7 +15,7 @@ namespace MbD { { public: void initialize() override; - void parseMBDyn(std::string line); + void parseMBDyn(std::string line) override; void readVelocity(std::vector& args); void readOmega(std::vector& args); diff --git a/OndselSolver/MBDynStructural.h b/OndselSolver/MBDynStructural.h index d38b3d3..542cd5b 100644 --- a/OndselSolver/MBDynStructural.h +++ b/OndselSolver/MBDynStructural.h @@ -14,7 +14,7 @@ namespace MbD { { public: MBDynStructural(); - void parseMBDyn(std::string line); + void parseMBDyn(std::string line) override; void readVelocity(std::vector& args); void readOmega(std::vector& args); void createASMT() override; From 7083905e59849b2bcc6e7c80ae229d85ad5896c9 Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Thu, 30 Nov 2023 22:25:48 -0700 Subject: [PATCH 03/12] MacOS edits 2 --- OndselSolver/ASMTAssembly.cpp | 2 +- OndselSolver/ASMTAssembly.h | 15 +++------------ OndselSolver/ASMTItem.cpp | 1 + OndselSolver/EulerParameters.ref.h | 6 ------ OndselSolver/OndselSolver.vcxproj | 1 - OndselSolver/OndselSolver.vcxproj.filters | 3 --- 6 files changed, 5 insertions(+), 23 deletions(-) delete mode 100644 OndselSolver/EulerParameters.ref.h diff --git a/OndselSolver/ASMTAssembly.cpp b/OndselSolver/ASMTAssembly.cpp index 0897438..830ed42 100644 --- a/OndselSolver/ASMTAssembly.cpp +++ b/OndselSolver/ASMTAssembly.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include "ASMTAssembly.h" #include "CREATE.h" @@ -40,7 +41,6 @@ #include "SystemSolver.h" #include "ASMTItemIJ.h" #include "ASMTKinematicIJ.h" -#include #include "ASMTAngleJoint.h" #include "ASMTConstantVelocityJoint.h" #include "ASMTCylSphJoint.h" diff --git a/OndselSolver/ASMTAssembly.h b/OndselSolver/ASMTAssembly.h index 28fff85..ed10d22 100644 --- a/OndselSolver/ASMTAssembly.h +++ b/OndselSolver/ASMTAssembly.h @@ -10,31 +10,22 @@ #include #include "ASMTSpatialContainer.h" -#include "FullColumn.h" -#include "FullMatrix.h" -#include "MBDynSystem.h" -#include "ASMTTime.h" -#include "ASMTConstantGravity.h" -#include "ASMTSimulationParameters.h" -#include "ASMTAnimationParameters.h" namespace MbD { - class ASMTRefPoint; - class ASMTRefCurve; - class ASMTRefSurface; class ASMTPart; class ASMTKinematicIJ; class ASMTConstraintSet; class ASMTForceTorque; class ASMTConstantGravity; - // class ASMTSimulationParameters; - // class ASMTAnimationParameters; + class ASMTSimulationParameters; + class ASMTAnimationParameters; class ASMTJoint; class ASMTMotion; class Units; class ASMTTime; class SystemSolver; class ASMTItemIJ; + class MBDynSystem; class ASMTAssembly : public ASMTSpatialContainer { diff --git a/OndselSolver/ASMTItem.cpp b/OndselSolver/ASMTItem.cpp index d41d546..9e40991 100644 --- a/OndselSolver/ASMTItem.cpp +++ b/OndselSolver/ASMTItem.cpp @@ -184,6 +184,7 @@ std::shared_ptr MbD::ASMTItem::sptrConstant(double value) void MbD::ASMTItem::storeOnLevel(std::ofstream& os, int level) { + noop(); assert(false); } diff --git a/OndselSolver/EulerParameters.ref.h b/OndselSolver/EulerParameters.ref.h deleted file mode 100644 index ca14f54..0000000 --- a/OndselSolver/EulerParameters.ref.h +++ /dev/null @@ -1,6 +0,0 @@ -#pragma once - -namespace MbD { - template - class EulerParameters; -} diff --git a/OndselSolver/OndselSolver.vcxproj b/OndselSolver/OndselSolver.vcxproj index 89ac539..24cfefc 100644 --- a/OndselSolver/OndselSolver.vcxproj +++ b/OndselSolver/OndselSolver.vcxproj @@ -539,7 +539,6 @@ - diff --git a/OndselSolver/OndselSolver.vcxproj.filters b/OndselSolver/OndselSolver.vcxproj.filters index 271f825..014e704 100644 --- a/OndselSolver/OndselSolver.vcxproj.filters +++ b/OndselSolver/OndselSolver.vcxproj.filters @@ -1814,9 +1814,6 @@ Header Files - - Header Files - Header Files From 9f2f1c4fd99912561ed9002bd1cf703f1917a7f3 Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Thu, 30 Nov 2023 22:42:46 -0700 Subject: [PATCH 04/12] Edits for MacOS 3 --- OndselSolver/ASMTAssembly.h | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/OndselSolver/ASMTAssembly.h b/OndselSolver/ASMTAssembly.h index ed10d22..7f24a76 100644 --- a/OndselSolver/ASMTAssembly.h +++ b/OndselSolver/ASMTAssembly.h @@ -10,19 +10,20 @@ #include #include "ASMTSpatialContainer.h" +//Required for initialization +#include "ASMTConstantGravity.h" +#include "ASMTSimulationParameters.h" +#include "ASMTAnimationParameters.h" +#include "ASMTTime.h" +#include "Units.h" namespace MbD { class ASMTPart; class ASMTKinematicIJ; class ASMTConstraintSet; class ASMTForceTorque; - class ASMTConstantGravity; - class ASMTSimulationParameters; - class ASMTAnimationParameters; class ASMTJoint; class ASMTMotion; - class Units; - class ASMTTime; class SystemSolver; class ASMTItemIJ; class MBDynSystem; From cafc91c71d737394853c21fa074ffbb608c8a21a Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Thu, 30 Nov 2023 23:12:17 -0700 Subject: [PATCH 05/12] EXPORT added --- OndselSolver/ASMTAssembly.h | 2 +- OndselSolver/ASMTInLineJoint.h | 2 +- OndselSolver/ASMTInPlaneJoint.h | 2 +- OndselSolver/ASMTItem.cpp | 1 + OndselSolver/ASMTKinematicIJ.h | 2 +- 5 files changed, 5 insertions(+), 4 deletions(-) diff --git a/OndselSolver/ASMTAssembly.h b/OndselSolver/ASMTAssembly.h index 7f24a76..8364634 100644 --- a/OndselSolver/ASMTAssembly.h +++ b/OndselSolver/ASMTAssembly.h @@ -28,7 +28,7 @@ namespace MbD { class ASMTItemIJ; class MBDynSystem; - class ASMTAssembly : public ASMTSpatialContainer + class EXPORT ASMTAssembly : public ASMTSpatialContainer { // public: diff --git a/OndselSolver/ASMTInLineJoint.h b/OndselSolver/ASMTInLineJoint.h index bac3829..371084f 100644 --- a/OndselSolver/ASMTInLineJoint.h +++ b/OndselSolver/ASMTInLineJoint.h @@ -11,7 +11,7 @@ #include "ASMTJoint.h" namespace MbD { - class ASMTInLineJoint : public ASMTJoint + class EXPORT ASMTInLineJoint : public ASMTJoint { // public: diff --git a/OndselSolver/ASMTInPlaneJoint.h b/OndselSolver/ASMTInPlaneJoint.h index d6bdf3d..445dcc8 100644 --- a/OndselSolver/ASMTInPlaneJoint.h +++ b/OndselSolver/ASMTInPlaneJoint.h @@ -11,7 +11,7 @@ #include "ASMTJoint.h" namespace MbD { - class ASMTInPlaneJoint : public ASMTJoint + class EXPORT ASMTInPlaneJoint : public ASMTJoint { // public: diff --git a/OndselSolver/ASMTItem.cpp b/OndselSolver/ASMTItem.cpp index 9e40991..4fd9c51 100644 --- a/OndselSolver/ASMTItem.cpp +++ b/OndselSolver/ASMTItem.cpp @@ -151,6 +151,7 @@ void MbD::ASMTItem::deleteMbD() void MbD::ASMTItem::createMbD(std::shared_ptr mbdSys, std::shared_ptr mbdUnits) { + noop(); assert(false); } diff --git a/OndselSolver/ASMTKinematicIJ.h b/OndselSolver/ASMTKinematicIJ.h index d6871e7..1948e57 100644 --- a/OndselSolver/ASMTKinematicIJ.h +++ b/OndselSolver/ASMTKinematicIJ.h @@ -11,7 +11,7 @@ #include "ASMTItemIJ.h" namespace MbD { - class ASMTKinematicIJ : public ASMTItemIJ + class EXPORT ASMTKinematicIJ : public ASMTItemIJ { // public: From 1e82769a8f3f067c41a8010abc5667a81294e417 Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Fri, 1 Dec 2023 00:02:17 -0700 Subject: [PATCH 06/12] storeOnTimeSeries --- OndselSolver/ASMTRevoluteJoint.cpp | 12 ++++++------ OndselSolver/ASMTRevoluteJoint.h | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/OndselSolver/ASMTRevoluteJoint.cpp b/OndselSolver/ASMTRevoluteJoint.cpp index 6193b49..4c3171b 100644 --- a/OndselSolver/ASMTRevoluteJoint.cpp +++ b/OndselSolver/ASMTRevoluteJoint.cpp @@ -16,9 +16,9 @@ std::shared_ptr MbD::ASMTRevoluteJoint::mbdClassNew() { return CREATE::With(); } - -void MbD::ASMTRevoluteJoint::storeOnTimeSeries(std::ofstream& os) -{ - os << "RevoluteJointSeries\t" << fullName("") << std::endl; - ASMTItemIJ::storeOnTimeSeries(os); -} +// +//void MbD::ASMTRevoluteJoint::storeOnTimeSeries(std::ofstream& os) +//{ +// os << "RevoluteJointSeries\t" << fullName("") << std::endl; +// ASMTItemIJ::storeOnTimeSeries(os); +//} diff --git a/OndselSolver/ASMTRevoluteJoint.h b/OndselSolver/ASMTRevoluteJoint.h index b506dfd..71cc274 100644 --- a/OndselSolver/ASMTRevoluteJoint.h +++ b/OndselSolver/ASMTRevoluteJoint.h @@ -16,7 +16,7 @@ namespace MbD { // public: std::shared_ptr mbdClassNew() override; - void storeOnTimeSeries(std::ofstream& os) override; + //void storeOnTimeSeries(std::ofstream& os) override; }; } From c11cbc2c72e63614115059bc01b09bec177b9a29 Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Fri, 1 Dec 2023 00:20:50 -0700 Subject: [PATCH 07/12] #include "ASMTKinematicIJ.h" --- OndselSolver/ASMTAssembly.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/OndselSolver/ASMTAssembly.cpp b/OndselSolver/ASMTAssembly.cpp index 830ed42..4e72a68 100644 --- a/OndselSolver/ASMTAssembly.cpp +++ b/OndselSolver/ASMTAssembly.cpp @@ -40,7 +40,6 @@ #include "ASMTTime.h" #include "SystemSolver.h" #include "ASMTItemIJ.h" -#include "ASMTKinematicIJ.h" #include "ASMTAngleJoint.h" #include "ASMTConstantVelocityJoint.h" #include "ASMTCylSphJoint.h" @@ -56,6 +55,7 @@ #include "ASMTRackPinionJoint.h" #include "ASMTScrewJoint.h" #include "SimulationStoppingError.h" +#include "ASMTKinematicIJ.h" using namespace MbD; From 2a3e026ee0a2fa7d6a9876f81862e3c767d7e90b Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Fri, 1 Dec 2023 00:36:06 -0700 Subject: [PATCH 08/12] move includes --- OndselSolver/ASMTAssembly.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/OndselSolver/ASMTAssembly.cpp b/OndselSolver/ASMTAssembly.cpp index 4e72a68..a2f9669 100644 --- a/OndselSolver/ASMTAssembly.cpp +++ b/OndselSolver/ASMTAssembly.cpp @@ -26,7 +26,6 @@ #include "ASMTFixedJoint.h" #include "ASMTGeneralMotion.h" #include "ASMTUniversalJoint.h" -#include "ExternalSystem.h" #include "ASMTPointInPlaneJoint.h" #include "ASMTPrincipalMassMarker.h" #include "ASMTForceTorque.h" @@ -34,11 +33,7 @@ #include "ASMTSimulationParameters.h" #include "ASMTAnimationParameters.h" #include "Part.h" -#include "ASMTRefPoint.h" -#include "ASMTRefCurve.h" -#include "ASMTRefSurface.h" #include "ASMTTime.h" -#include "SystemSolver.h" #include "ASMTItemIJ.h" #include "ASMTAngleJoint.h" #include "ASMTConstantVelocityJoint.h" @@ -56,6 +51,11 @@ #include "ASMTScrewJoint.h" #include "SimulationStoppingError.h" #include "ASMTKinematicIJ.h" +#include "ASMTRefPoint.h" +#include "ASMTRefCurve.h" +#include "ASMTRefSurface.h" +#include "ExternalSystem.h" +#include "SystemSolver.h" using namespace MbD; From a2871455f418fb75ddff1d95db75ad2c69ce5173 Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Fri, 1 Dec 2023 00:53:36 -0700 Subject: [PATCH 09/12] ASMTPointInPlaneJoint::storeOnLevel --- OndselSolver/ASMTPointInPlaneJoint.cpp | 5 +++++ OndselSolver/ASMTPointInPlaneJoint.h | 1 + 2 files changed, 6 insertions(+) diff --git a/OndselSolver/ASMTPointInPlaneJoint.cpp b/OndselSolver/ASMTPointInPlaneJoint.cpp index 733055a..3a195fc 100644 --- a/OndselSolver/ASMTPointInPlaneJoint.cpp +++ b/OndselSolver/ASMTPointInPlaneJoint.cpp @@ -23,3 +23,8 @@ void MbD::ASMTPointInPlaneJoint::storeOnTimeSeries(std::ofstream& os) ASMTItemIJ::storeOnTimeSeries(os); } +void MbD::ASMTPointInPlaneJoint::storeOnLevel(std::ofstream& os, int level) +{ + ASMTInPlaneJoint::storeOnLevel(os, level); +} + diff --git a/OndselSolver/ASMTPointInPlaneJoint.h b/OndselSolver/ASMTPointInPlaneJoint.h index df1e903..ad31541 100644 --- a/OndselSolver/ASMTPointInPlaneJoint.h +++ b/OndselSolver/ASMTPointInPlaneJoint.h @@ -17,6 +17,7 @@ namespace MbD { public: std::shared_ptr mbdClassNew() override; void storeOnTimeSeries(std::ofstream& os) override; + void storeOnLevel(std::ofstream& os, int level) override; }; } From 428e81c5cec63783114b8cd3c687f9816a2b3308 Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Fri, 1 Dec 2023 01:13:13 -0700 Subject: [PATCH 10/12] using namespace MbD; --- OndselSolver/ASMTAtPointJoint.cpp | 2 ++ OndselSolver/ASMTInLineJoint.cpp | 2 ++ OndselSolver/ASMTInPlaneJoint.cpp | 2 ++ 3 files changed, 6 insertions(+) diff --git a/OndselSolver/ASMTAtPointJoint.cpp b/OndselSolver/ASMTAtPointJoint.cpp index b43f8db..f0b3efe 100644 --- a/OndselSolver/ASMTAtPointJoint.cpp +++ b/OndselSolver/ASMTAtPointJoint.cpp @@ -1 +1,3 @@ #include "ASMTAtPointJoint.h" + +using namespace MbD; diff --git a/OndselSolver/ASMTInLineJoint.cpp b/OndselSolver/ASMTInLineJoint.cpp index 27e9382..5aee6f6 100644 --- a/OndselSolver/ASMTInLineJoint.cpp +++ b/OndselSolver/ASMTInLineJoint.cpp @@ -1,2 +1,4 @@ #include "ASMTInLineJoint.h" + +using namespace MbD; diff --git a/OndselSolver/ASMTInPlaneJoint.cpp b/OndselSolver/ASMTInPlaneJoint.cpp index 1b372da..b15371d 100644 --- a/OndselSolver/ASMTInPlaneJoint.cpp +++ b/OndselSolver/ASMTInPlaneJoint.cpp @@ -2,6 +2,8 @@ #include "ASMTInPlaneJoint.h" #include "InPlaneJoint.h" +using namespace MbD; + void MbD::ASMTInPlaneJoint::parseASMT(std::vector& lines) { ASMTJoint::parseASMT(lines); From 173b94323a54af1bb60adb5eee3e1548bfba5286 Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Fri, 1 Dec 2023 01:39:50 -0700 Subject: [PATCH 11/12] MBDynGravity --- OndselSolver/ASMTPointInPlaneJoint.cpp | 5 ----- OndselSolver/ASMTPointInPlaneJoint.h | 1 - OndselSolver/MBDynGravity.cpp | 1 + 3 files changed, 1 insertion(+), 6 deletions(-) diff --git a/OndselSolver/ASMTPointInPlaneJoint.cpp b/OndselSolver/ASMTPointInPlaneJoint.cpp index 3a195fc..733055a 100644 --- a/OndselSolver/ASMTPointInPlaneJoint.cpp +++ b/OndselSolver/ASMTPointInPlaneJoint.cpp @@ -23,8 +23,3 @@ void MbD::ASMTPointInPlaneJoint::storeOnTimeSeries(std::ofstream& os) ASMTItemIJ::storeOnTimeSeries(os); } -void MbD::ASMTPointInPlaneJoint::storeOnLevel(std::ofstream& os, int level) -{ - ASMTInPlaneJoint::storeOnLevel(os, level); -} - diff --git a/OndselSolver/ASMTPointInPlaneJoint.h b/OndselSolver/ASMTPointInPlaneJoint.h index ad31541..df1e903 100644 --- a/OndselSolver/ASMTPointInPlaneJoint.h +++ b/OndselSolver/ASMTPointInPlaneJoint.h @@ -17,7 +17,6 @@ namespace MbD { public: std::shared_ptr mbdClassNew() override; void storeOnTimeSeries(std::ofstream& os) override; - void storeOnLevel(std::ofstream& os, int level) override; }; } diff --git a/OndselSolver/MBDynGravity.cpp b/OndselSolver/MBDynGravity.cpp index 02b0f76..a81a9a6 100644 --- a/OndselSolver/MBDynGravity.cpp +++ b/OndselSolver/MBDynGravity.cpp @@ -48,6 +48,7 @@ void MbD::MBDynGravity::parseMBDyn(std::string line) void MbD::MBDynGravity::readFunction(std::vector& args) { assert(false); + noop(); } void MbD::MBDynGravity::createASMT() From 9266d9cdf93b067a2aec6f4108f4e0852ba2fd78 Mon Sep 17 00:00:00 2001 From: Aik-Siong Koh Date: Fri, 1 Dec 2023 11:25:54 -0700 Subject: [PATCH 12/12] CMakeLists.txt --- CMakeLists.txt | 30 +++++++++++++++++++++++ OndselSolver/OndselSolver.vcxproj | 3 +++ OndselSolver/OndselSolver.vcxproj.filters | 3 +++ 3 files changed, 36 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 12e4c98..589d6ff 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,36 +46,50 @@ set(ONDSELSOLVER_SRC OndselSolver/ArcSine.cpp OndselSolver/ArcTan.cpp OndselSolver/ArcTan2.cpp + OndselSolver/ASMTAngleJoint.cpp OndselSolver/ASMTAnimationParameters.cpp OndselSolver/ASMTAssembly.cpp + OndselSolver/ASMTCompoundJoint.cpp OndselSolver/ASMTConstantGravity.cpp + OndselSolver/ASMTConstantVelocityJoint.cpp OndselSolver/ASMTConstraintSet.cpp OndselSolver/ASMTCylindricalJoint.cpp + OndselSolver/ASMTCylSphJoint.cpp OndselSolver/ASMTExtrusion.cpp OndselSolver/ASMTFixedJoint.cpp OndselSolver/ASMTForceTorque.cpp + OndselSolver/ASMTGearJoint.cpp OndselSolver/ASMTGeneralMotion.cpp + OndselSolver/ASMTInPlaneJoint.cpp OndselSolver/ASMTItem.cpp OndselSolver/ASMTItemIJ.cpp OndselSolver/ASMTJoint.cpp OndselSolver/ASMTKinematicIJ.cpp + OndselSolver/ASMTLineInPlaneJoint.cpp OndselSolver/ASMTMarker.cpp OndselSolver/ASMTMotion.cpp OndselSolver/ASMTNoRotationJoint.cpp + OndselSolver/ASMTParallelAxesJoint.cpp OndselSolver/ASMTPart.cpp + OndselSolver/ASMTPerpendicularJoint.cpp + OndselSolver/ASMTPlanarJoint.cpp OndselSolver/ASMTPointInLineJoint.cpp OndselSolver/ASMTPointInPlaneJoint.cpp OndselSolver/ASMTPrincipalMassMarker.cpp + OndselSolver/ASMTRackPinionJoint.cpp OndselSolver/ASMTRefCurve.cpp OndselSolver/ASMTRefItem.cpp OndselSolver/ASMTRefPoint.cpp OndselSolver/ASMTRefSurface.cpp + OndselSolver/ASMTRevCylJoint.cpp OndselSolver/ASMTRevoluteJoint.cpp OndselSolver/ASMTRotationalMotion.cpp + OndselSolver/ASMTScrewJoint.cpp OndselSolver/ASMTSimulationParameters.cpp OndselSolver/ASMTSpatialContainer.cpp OndselSolver/ASMTSpatialItem.cpp OndselSolver/ASMTSphericalJoint.cpp + OndselSolver/ASMTSphSphJoint.cpp OndselSolver/ASMTTime.cpp OndselSolver/ASMTTranslationalJoint.cpp OndselSolver/ASMTTranslationalMotion.cpp @@ -223,6 +237,7 @@ set(ONDSELSOLVER_SRC OndselSolver/MBDynData.cpp OndselSolver/MBDynDrive.cpp OndselSolver/MBDynElement.cpp + OndselSolver/MBDynGravity.cpp OndselSolver/MBDynInitialValue.cpp OndselSolver/MBDynItem.cpp OndselSolver/MBDynJoint.cpp @@ -338,36 +353,50 @@ set(ONDSELSOLVER_HEADERS OndselSolver/ArcSine.h OndselSolver/ArcTan.h OndselSolver/ArcTan2.h + OndselSolver/ASMTAngleJoint.h OndselSolver/ASMTAnimationParameters.h OndselSolver/ASMTAssembly.h + OndselSolver/ASMTCylSphJoint.h + OndselSolver/ASMTCompoundJoint.h OndselSolver/ASMTConstantGravity.h + OndselSolver/ASMTConstantVelocityJoint.h OndselSolver/ASMTConstraintSet.h OndselSolver/ASMTCylindricalJoint.h OndselSolver/ASMTExtrusion.h OndselSolver/ASMTFixedJoint.h OndselSolver/ASMTForceTorque.h + OndselSolver/ASMTGearJoint.h OndselSolver/ASMTGeneralMotion.h + OndselSolver/ASMTInPlaneJoint.h OndselSolver/ASMTItem.h OndselSolver/ASMTItemIJ.h OndselSolver/ASMTJoint.h OndselSolver/ASMTKinematicIJ.h + OndselSolver/ASMTLineInPlaneJoint.h OndselSolver/ASMTMarker.h OndselSolver/ASMTMotion.h OndselSolver/ASMTNoRotationJoint.h + OndselSolver/ASMTParallelAxesJoint.h OndselSolver/ASMTPart.h + OndselSolver/ASMTPerpendicularJoint.h + OndselSolver/ASMTPlanarJoint.h OndselSolver/ASMTPointInLineJoint.h OndselSolver/ASMTPointInPlaneJoint.h OndselSolver/ASMTPrincipalMassMarker.h + OndselSolver/ASMTRackPinionJoint.h OndselSolver/ASMTRefCurve.h OndselSolver/ASMTRefItem.h OndselSolver/ASMTRefPoint.h OndselSolver/ASMTRefSurface.h + OndselSolver/ASMTRevCylJoint.h OndselSolver/ASMTRevoluteJoint.h OndselSolver/ASMTRotationalMotion.h + OndselSolver/ASMTScrewJoint.h OndselSolver/ASMTSimulationParameters.h OndselSolver/ASMTSpatialContainer.h OndselSolver/ASMTSpatialItem.h OndselSolver/ASMTSphericalJoint.h + OndselSolver/ASMTSphSphJoint.h OndselSolver/ASMTTime.h OndselSolver/ASMTTranslationalJoint.h OndselSolver/ASMTTranslationalMotion.h @@ -516,6 +545,7 @@ set(ONDSELSOLVER_HEADERS OndselSolver/MBDynControlData.h OndselSolver/MBDynData.h OndselSolver/MBDynElement.h + OndselSolver/MBDynGravity.h OndselSolver/MBDynInitialValue.h OndselSolver/MBDynItem.h OndselSolver/MBDynJoint.h diff --git a/OndselSolver/OndselSolver.vcxproj b/OndselSolver/OndselSolver.vcxproj index 24cfefc..a869c2f 100644 --- a/OndselSolver/OndselSolver.vcxproj +++ b/OndselSolver/OndselSolver.vcxproj @@ -797,6 +797,9 @@ + + + diff --git a/OndselSolver/OndselSolver.vcxproj.filters b/OndselSolver/OndselSolver.vcxproj.filters index 014e704..1e5a4bb 100644 --- a/OndselSolver/OndselSolver.vcxproj.filters +++ b/OndselSolver/OndselSolver.vcxproj.filters @@ -1912,4 +1912,7 @@ + + + \ No newline at end of file