diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 5960aa0..1333e80 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -45,10 +45,11 @@ void DRMSimulationModel::runUnitTests() Settings settings; settings.Dtini = 0.1; settings.xi = 0.9969; - settings.maxDRMIterations = 0; - formFinder.mSettings.totalResidualForcesNormThreshold = 1; + settings.totalResidualForcesNormThreshold = 1; settings.shouldDraw = false; settings.beVerbose = true; + // settings.debugModeStep = 1000; + // settings.shouldDraw = true; settings.shouldCreatePlots = true; SimulationResults simpleBeam_simulationResults = formFinder.executeSimulation(std::make_shared(beamSimulationJob), settings); @@ -104,10 +105,11 @@ void DRMSimulationModel::runUnitTests() shortSpanGridshellSimulationJob.pMesh->setBeamMaterial(0.3, 200 * 1e9); assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); shortSpanGridshellSimulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026}); + DRMSimulationModel formFinder2; SimulationResults shortSpanGridshellSimulationResults - = formFinder.executeSimulation(std::make_shared( - shortSpanGridshellSimulationJob), - settings); + = formFinder2.executeSimulation(std::make_shared( + shortSpanGridshellSimulationJob), + settings); shortSpanGridshellSimulationResults.save(); const std::string groundOfTruthBinaryFilename @@ -183,10 +185,11 @@ void DRMSimulationModel::runUnitTests() } longSpanGridshell_simulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026}); + DRMSimulationModel formFinder3; SimulationResults longSpanGridshell_simulationResults - = formFinder.executeSimulation(std::make_shared( - longSpanGridshell_simulationJob), - settings); + = formFinder3.executeSimulation(std::make_shared( + longSpanGridshell_simulationJob), + settings); longSpanGridshell_simulationResults.save(); const std::string longSpan_groundOfTruthBinaryFilename @@ -217,22 +220,14 @@ void DRMSimulationModel::runUnitTests() // polyscope::show(); } - void DRMSimulationModel::reset(const std::shared_ptr &pJob) { - mCurrentSimulationStep = 0; - history.clear(); - history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel(); - plotYValues.clear(); - plotHandle.reset(); - checkedForMaximumMoment = false; - externalMomentsNorm = 0; - Dt = mSettings.Dtini; - numOfDampings = 0; - shouldTemporarilyDampForces = false; - externalLoadStep = 1; - isVertexConstrained.clear(); - minTotalResidualForcesNorm = std::numeric_limits::max(); + //#ifdef USE_ENSMALLEN + // this->pJob = pJob; + //#endif + pMesh.reset(); + pMesh = std::make_unique(*pJob->pMesh); + vcg::tri::UpdateBounding::Box(*pMesh); constrainedVertices.clear(); constrainedVertices = pJob->constrainedVertices; @@ -247,6 +242,26 @@ void DRMSimulationModel::reset(const std::shared_ptr &pJob) for (auto fixedVertex : constrainedVertices) { isVertexConstrained[fixedVertex.first] = true; } +} + +void DRMSimulationModel::reset(const std::shared_ptr &pJob, const Settings &settings) +{ + mSettings = settings; + mCurrentSimulationStep = 0; + history.clear(); + history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel(); + plotYValues.clear(); + plotHandle.reset(); + checkedForMaximumMoment = false; + externalMomentsNorm = 0; + Dt = settings.Dtini; + numOfDampings = 0; + shouldTemporarilyDampForces = false; + externalLoadStep = 1; + isVertexConstrained.clear(); + minTotalResidualForcesNorm = std::numeric_limits::max(); + + reset(pJob); #ifdef POLYSCOPE_DEFINED if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) { @@ -257,14 +272,14 @@ void DRMSimulationModel::reset(const std::shared_ptr &pJob) polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel + "_" + pMesh->getLabel(), pMesh->getEigenVertices(), pMesh->getEigenEdges()) - ->setRadius(pMesh->elements[0].dimensions.getDrawingRadius()); + ->setRadius(0.002); // registerWorldAxes(); } #endif - // if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) { - // shouldApplyInitialDistortion = true; - // } + if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) { + shouldApplyInitialDistortion = true; + } updateElementalFrames(); } @@ -1414,7 +1429,6 @@ void DRMSimulationModel::updateResidualForces() void DRMSimulationModel::computeRigidSupports() { - assert(pMesh != nullptr); isRigidSupport.clear(); isRigidSupport.resize(pMesh->VN(), false); for (const VertexType &v : pMesh->vert) { @@ -1625,6 +1639,7 @@ void DRMSimulationModel::updateNodalVelocities() for (VertexType &v : pMesh->vert) { const VertexIndex vi = pMesh->getIndex(v); Node &node = pMesh->nodes[v]; + node.previousVelocity = node.velocity; if (mSettings.viscousDampingFactor.has_value()) { const Vector6d massOverDt = node.mass_6d / Dt; // const Vector6d visciousDampingFactor(viscuousDampingConstant / 2); @@ -1777,11 +1792,17 @@ void DRMSimulationModel::updateNodeNormal( void DRMSimulationModel::applyDisplacements( const std::unordered_map> &fixedVertices) { - for (VertexType &v : pMesh->vert) { - updateNodePosition(v, fixedVertices); - updateNodeNormal(v, fixedVertices); - updateNodeNr(v); - } + std::for_each( +#ifdef ENABLE_PARALLEL_DRM + std::execution::par_unseq, +#endif + pMesh->vert.begin(), + pMesh->vert.end(), + [&](auto &v) { + updateNodePosition(v, fixedVertices); + updateNodeNormal(v, fixedVertices); + updateNodeNr(v); + }); updateElementalFrames(); if (mSettings.shouldDraw) { pMesh->updateEigenEdgeAndVertices(); @@ -1847,8 +1868,12 @@ void DRMSimulationModel::applyForcedNormals( // NOTE: Is this correct? Should the kinetic energy be computed like that? void DRMSimulationModel::updateKineticEnergy() { + pMesh->pre_previousTotalKineticEnergy = pMesh->previousTotalKineticEnergy; + pMesh->pre_previousTotalTranslationalKineticEnergy + = pMesh->previousTotalTranslationalKineticEnergy; + pMesh->pre_previousTotalRotationalKineticEnergy = pMesh->previousTotalRotationalKineticEnergy; pMesh->previousTotalKineticEnergy = pMesh->currentTotalKineticEnergy; - pMesh->previousTranslationalKineticEnergy = pMesh->currentTotalTranslationalKineticEnergy; + pMesh->previousTotalTranslationalKineticEnergy = pMesh->currentTotalTranslationalKineticEnergy; pMesh->previousTotalRotationalKineticEnergy = pMesh->currentTotalRotationalKineticEnergy; pMesh->currentTotalKineticEnergy = 0; pMesh->currentTotalTranslationalKineticEnergy = 0; @@ -1889,7 +1914,14 @@ void DRMSimulationModel::resetVelocities() // // reset // // current to 0 or this? 0; + pMesh->nodes[v].previousVelocity = 0; } + pMesh->pre_previousTotalKineticEnergy = 0; + pMesh->pre_previousTotalTranslationalKineticEnergy = 0; + pMesh->pre_previousTotalRotationalKineticEnergy = 0; + pMesh->previousTotalKineticEnergy = 0; + pMesh->previousTotalTranslationalKineticEnergy = 0; + pMesh->previousTotalRotationalKineticEnergy = 0; updateKineticEnergy(); } @@ -1960,6 +1992,7 @@ void DRMSimulationModel::updatePositionsOnTheFly( for (VertexType &v : pMesh->vert) { Node &node = pMesh->nodes[v]; + node.previousVelocity = node.velocity; node.velocity = node.velocity + node.acceleration * Dt; } updateKineticEnergy(); @@ -2564,15 +2597,26 @@ void DRMSimulationModel::applyKineticDamping(const std::shared_ptrpre_previousTotalKineticEnergy; + // const double &KE2 = pMesh->previousTotalKineticEnergy; + // const double &KE3 = pMesh->currentTotalKineticEnergy; + const double &KE1 = pMesh->pre_previousTotalTranslationalKineticEnergy; + const double &KE2 = pMesh->previousTotalTranslationalKineticEnergy; + const double &KE3 = pMesh->currentTotalTranslationalKineticEnergy; + const double bitaDt = 0.5 * Dt * (3 + (-3 * KE1 + 4 * KE2 - KE3) / (KE1 - 2 * KE2 + KE3)); for (VertexType &v : pMesh->vert) { Node &node = pMesh->nodes[v]; - Vector6d stepDisplacement = node.velocity * 0.5 * Dt; + // Vector6d stepDisplacement = node.velocity * 0.5 * Dt; // if (shouldCapDisplacements // && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) { // stepDisplacement = stepDisplacement // * (*mSettings.displacementCap // / stepDisplacement.getTranslation().norm()); // } + + // Vector6d stepDisplacement = node.velocity * 0.5 * Dt; + Vector6d stepDisplacement = node.velocity * Dt + node.previousVelocity * bitaDt; + node.displacements = node.displacements - stepDisplacement; } applyDisplacements(constrainedVertices); @@ -2598,51 +2642,38 @@ void DRMSimulationModel::applyKineticDamping(const std::shared_ptr &pJob) { - assert(pMesh != nullptr && pMesh->VN() == pJob->pMesh->VN() && pMesh->EN() == pJob->pMesh->EN()); - reset(pJob); auto beginTime = std::chrono::high_resolution_clock::now(); updateNodalMasses(); - // constexpr bool useDRM = true; - //#ifdef USE_ENSMALLEN - // if (!useDRM) { - // setJob(pJob); - // // ens::L_BFGS optimizer(20); - // ens::SA optimizer; - // arma::mat x(pJob->pMesh->VN() * NumDoF, 1); - // optimizer.Optimize(*this, x); - // // getD - // } else { - //#endif - // std::unordered_map nodalExternalForces = pJob->nodalExternalForces; - // double totalExternalForcesNorm = 0; - // Vector6d sumOfExternalForces(0); - // for (auto &nodalForce : nodalExternalForces) { - // const double percentageOfExternalLoads = double(externalLoadStep) - // / mSettings.desiredGradualExternalLoadsSteps; - // nodalForce.second = nodalForce.second * percentageOfExternalLoads; - // totalExternalForcesNorm += nodalForce.second.norm(); - // // sumOfExternalForces = sumOfExternalForces + nodalForce.second; - // } updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices); - if (!pJob->nodalExternalForces.empty()) { + if (!pJob->nodalExternalForces.empty() + && !mSettings.totalResidualForcesNormThreshold.has_value()) { mSettings.totalResidualForcesNormThreshold = mSettings.totalExternalForcesNormPercentageTermination * pMesh->totalExternalForcesNorm; - } else { + } /*else { mSettings.totalResidualForcesNormThreshold = 1e-3; std::cout << "No forces setted default residual forces norm threshold" << std::endl; - } - if (mSettings.beVerbose) { - std::cout << "totalResidualForcesNormThreshold:" - << mSettings.totalResidualForcesNormThreshold << std::endl; - if (mSettings.averageResidualForcesCriterionThreshold.has_value()) { - std::cout << "average/extNorm threshold:" - << *mSettings.averageResidualForcesCriterionThreshold << std::endl; - } + }*/ + if (mSettings.beVerbose) { + // std::cout << "totalResidualForcesNormThreshold:" + // << mSettings.totalResidualForcesNormThreshold << std::endl; + if (mSettings.averageResidualForcesCriterionThreshold.has_value()) { + std::cout << "average/extNorm threshold:" + << *mSettings.averageResidualForcesCriterionThreshold << std::endl; } + } if (mSettings.beVerbose) { std::cout << "Executing simulation for mesh with:" << pMesh->VN() << " nodes and " @@ -2655,25 +2686,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrsave("./PatternOptimizationNonConv"); - // Dt = mSettings.Dtini; - // } - // if (mCurrentSimulationStep == 500 && shouldTemporarilyDampForces) { - // Dt = mSettings.Dtini; - // } - // while (true) { - updateNormalDerivatives(); - updateT1Derivatives(); - updateRDerivatives(); - updateT2Derivatives(); - updateT3Derivatives(); - const bool shouldBreak = mCurrentSimulationStep == 3935; - if (shouldBreak) { - int i = 0; - i++; - } + updateDerivatives(); updateResidualForcesOnTheFly(constrainedVertices); // TODO: write parallel function for updating positions @@ -2690,9 +2703,6 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrnodalForcedDisplacements); } - // if (!pJob->nodalForcedNormals.empty()) { - // applyForcedNormals(pJob->nodalForcedNormals); - // } updateElementalLengths(); mCurrentSimulationStep++; if (std::isnan(pMesh->currentTotalKineticEnergy)) { @@ -2790,8 +2800,8 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrtotalResidualForcesNorm); + // percentageOfConvergence.push_back(100 * mSettings.totalResidualForcesNormThreshold + // / pMesh->totalResidualForcesNorm); } if (mSettings.shouldCreatePlots && mSettings.debugModeStep.has_value() @@ -2880,15 +2890,19 @@ currentSimulationStep > maxDRMIterations*/ // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm // << std::endl; const bool fullfillsResidualForcesNormTerminationCriterion - = !mSettings.averageResidualForcesCriterionThreshold.has_value() + = mSettings.averageResidualForcesCriterionThreshold.has_value() && pMesh->totalResidualForcesNorm / pMesh->totalExternalForcesNorm < mSettings.totalExternalForcesNormPercentageTermination; const bool fullfillsAverageResidualForcesNormTerminationCriterion = mSettings.averageResidualForcesCriterionThreshold.has_value() && (pMesh->totalResidualForcesNorm / pMesh->VN()) / pMesh->totalExternalForcesNorm < mSettings.averageResidualForcesCriterionThreshold.value(); + const bool fullfillsResidualForcesNormThreshold + = mSettings.totalResidualForcesNormThreshold.has_value() + && pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold; if ((fullfillsAverageResidualForcesNormTerminationCriterion - || fullfillsResidualForcesNormTerminationCriterion) + || fullfillsResidualForcesNormTerminationCriterion + || fullfillsResidualForcesNormThreshold) && numOfDampings > 0 && (pJob->nodalForcedDisplacements.empty() || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) { @@ -2906,12 +2920,15 @@ currentSimulationStep > maxDRMIterations*/ break; } - const bool isKineticEnergyPeak = (mSettings.useTotalRotationalKineticEnergyForKineticDamping - && pMesh->previousTotalRotationalKineticEnergy - > pMesh->currentTotalRotationalKineticEnergy) - || (mSettings.useTranslationalKineticEnergyForKineticDamping - && pMesh->previousTranslationalKineticEnergy - > pMesh->currentTotalTranslationalKineticEnergy); + const bool isKineticEnergyPeak + = /* pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy + ||*/ + (mSettings.useTotalRotationalKineticEnergyForKineticDamping + && pMesh->previousTotalRotationalKineticEnergy + > pMesh->currentTotalRotationalKineticEnergy) + || (mSettings.useTranslationalKineticEnergyForKineticDamping + && pMesh->previousTotalTranslationalKineticEnergy + > pMesh->currentTotalTranslationalKineticEnergy); if (isKineticEnergyPeak) { const bool fullfillsKineticEnergyTerminationCriterion = mSettings.totalTranslationalKineticEnergyThreshold.has_value() @@ -2964,13 +2981,13 @@ currentSimulationStep > maxDRMIterations*/ // } } - if (mSettings.useTranslationalKineticEnergyForKineticDamping) { - applyKineticDamping(pJob); - if (mSettings.shouldCreatePlots) { - history.redMarks.push_back(mCurrentSimulationStep); - } - } + // if (mSettings.useTranslationalKineticEnergyForKineticDamping) { + applyKineticDamping(pJob); Dt *= mSettings.xi; + if (mSettings.shouldCreatePlots) { + history.redMarks.push_back(mCurrentSimulationStep); + } + // } // if (mSettings.isDebugMode) { // std::cout << Dt << std::endl; // } @@ -3002,11 +3019,9 @@ currentSimulationStep > maxDRMIterations*/ void DRMSimulationModel::setStructure(const std::shared_ptr &pMesh) { - this->pMesh.reset(); - this->pMesh = std::make_unique(*pMesh); - vcg::tri::UpdateBounding::Box(*pMesh); - // assert(false); - // std::terminate(); + std::cout << "This function is currently not implemented" << std::endl; + assert(false); + std::terminate(); } SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr &pJob, @@ -3014,8 +3029,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrpMesh); + reset(pJob, settings); assert(pJob->pMesh != nullptr); if (!solutionGuess.displacements.empty()) { diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 0267a6f..b7e4f2a 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -47,8 +47,8 @@ public: int gradualForcedDisplacementSteps{50}; // int desiredGradualExternalLoadsSteps{1}; double gamma{0.8}; - double totalResidualForcesNormThreshold{1e-20}; - double totalExternalForcesNormPercentageTermination{1e-3}; + std::optional totalResidualForcesNormThreshold; + double totalExternalForcesNormPercentageTermination{1e-5}; std::optional maxDRMIterations; std::optional debugModeStep; std::optional totalTranslationalKineticEnergyThreshold; @@ -104,10 +104,11 @@ private: SimulationHistory history; // Eigen::Tensor theta3Derivatives; // std::unordered_map theta3Derivatives; - bool shouldApplyInitialDistortion = false; + bool shouldApplyInitialDistortion{false}; //#ifdef USE_ENSMALLEN // std::shared_ptr pJob; //#endif + void reset(const std::shared_ptr &pJob, const Settings &settings); void updateNodalInternalForces( const std::unordered_map> &fixedVertices); void updateNodalExternalForces( @@ -259,6 +260,8 @@ private: std::vector> computeInternalForces( const std::unordered_map> &fixedVertices); + void updateDerivatives(); + public: DRMSimulationModel(); SimulationResults executeSimulation(const std::shared_ptr &pJob, diff --git a/simulationmesh.cpp b/simulationmesh.cpp index 6be0423..a3ba074 100755 --- a/simulationmesh.cpp +++ b/simulationmesh.cpp @@ -268,12 +268,12 @@ void SimulationMesh::unregister() const } #endif -void SimulationMesh::setBeamCrossSection( - const CrossSectionType &beamDimensions) { - for (size_t ei = 0; ei < EN(); ei++) { - elements[ei].dimensions = beamDimensions; - elements[ei].updateRigidity(); - } +void SimulationMesh::setBeamCrossSection(const CrossSectionType &beamDimensions) +{ + for (size_t ei = 0; ei < EN(); ei++) { + elements[ei].dimensions = beamDimensions; + elements[ei].updateRigidity(); + } } void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) { diff --git a/simulationmesh.hpp b/simulationmesh.hpp index 5ac43cb..21374e8 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -8,8 +8,8 @@ //extern bool drawGlobal; struct Element; struct Node; -using CrossSectionType = RectangularBeamDimensions; -//using CrossSectionType = CylindricalBeamDimensions; +//using CrossSectionType = RectangularBeamDimensions; +using CrossSectionType = CylindricalBeamDimensions; class SimulationMesh : public VCGEdgeMesh { private: @@ -38,8 +38,11 @@ public: std::vector getBeamDimensions(); std::vector getBeamMaterial(); + double pre_previousTotalKineticEnergy{0}; + double pre_previousTotalTranslationalKineticEnergy{0}; + double pre_previousTotalRotationalKineticEnergy{0}; double previousTotalKineticEnergy{0}; - double previousTranslationalKineticEnergy{0}; + double previousTotalTranslationalKineticEnergy{0}; double previousTotalRotationalKineticEnergy{0}; double previousTotalResidualForcesNorm{0}; double currentTotalKineticEnergy{0}; @@ -68,73 +71,76 @@ public: #endif }; -struct Element { +struct Element +{ + CrossSectionType dimensions; + ElementMaterial material; - CrossSectionType dimensions; - ElementMaterial material; + void computeMaterialProperties(const ElementMaterial &material); + // void computeDimensionsProperties(const RectangularBeamDimensions &dimensions); + // void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions); + void setDimensions(const CrossSectionType &dimensions); + void setMaterial(const ElementMaterial &material); + double getMass(const double &matrialDensity); - void computeMaterialProperties(const ElementMaterial &material); - // void computeDimensionsProperties(const RectangularBeamDimensions &dimensions); - // void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions); - void setDimensions(const CrossSectionType &dimensions); - void setMaterial(const ElementMaterial &material); - double getMass(const double &matrialDensity); + struct LocalFrame + { + VectorType t1; + VectorType t2; + VectorType t3; + }; - struct LocalFrame { - VectorType t1; - VectorType t2; - VectorType t3; - }; + EdgeIndex ei; + double length{0}; + double initialLength; + LocalFrame frame; - EdgeIndex ei; - double length{0}; - double initialLength; - LocalFrame frame; + struct Rigidity + { + double axial; + double torsional; + double firstBending; + double secondBending; + std::string toString() const + { + return std::string("Rigidity:") + std::string("\nAxial=") + std::to_string(axial) + + std::string("\nTorsional=") + std::to_string(torsional) + + std::string("\nFirstBending=") + std::to_string(firstBending) + + std::string("\nSecondBending=") + std::to_string(secondBending); + } + }; + Rigidity rigidity; + void updateRigidity(); - struct Rigidity { - double axial; - double torsional; - double firstBending; - double secondBending; - std::string toString() const { - return std::string("Rigidity:") + std::string("\nAxial=") + - std::to_string(axial) + std::string("\nTorsional=") + - std::to_string(torsional) + std::string("\nFirstBending=") + - std::to_string(firstBending) + std::string("\nSecondBending=") + - std::to_string(secondBending); - } - }; - Rigidity rigidity; - void updateRigidity(); + VectorType f1_j; + VectorType f1_jplus1; + VectorType f2_j; + VectorType f2_jplus1; + VectorType f3_j; + VectorType f3_jplus1; + double cosRotationAngle_j; + double cosRotationAngle_jplus1; + double sinRotationAngle_j; + double sinRotationAngle_jplus1; + std::vector> derivativeT1; + std::vector> derivativeT2; + std::vector> derivativeT3; + std::vector derivativeT1_j; + std::vector derivativeT1_jplus1; + std::vector derivativeT2_j; + std::vector derivativeT2_jplus1; + std::vector derivativeT3_j; + std::vector derivativeT3_jplus1; + std::vector derivativeR_j; + std::vector derivativeR_jplus1; + struct RotationalDisplacements + { + double theta1{0}, theta2{0}, theta3{0}; + }; + RotationalDisplacements rotationalDisplacements_j; + RotationalDisplacements rotationalDisplacements_jplus1; - VectorType f1_j; - VectorType f1_jplus1; - VectorType f2_j; - VectorType f2_jplus1; - VectorType f3_j; - VectorType f3_jplus1; - double cosRotationAngle_j; - double cosRotationAngle_jplus1; - double sinRotationAngle_j; - double sinRotationAngle_jplus1; - std::vector> derivativeT1; - std::vector> derivativeT2; - std::vector> derivativeT3; - std::vector derivativeT1_j; - std::vector derivativeT1_jplus1; - std::vector derivativeT2_j; - std::vector derivativeT2_jplus1; - std::vector derivativeT3_j; - std::vector derivativeT3_jplus1; - std::vector derivativeR_j; - std::vector derivativeR_jplus1; - struct RotationalDisplacements { - double theta1{0}, theta2{0}, theta3{0}; - }; - RotationalDisplacements rotationalDisplacements_j; - RotationalDisplacements rotationalDisplacements_jplus1; - - static void computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A); + static void computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A); }; struct Node { @@ -164,6 +170,7 @@ struct Node { CoordType initialNormal; Vector6d acceleration{0}; Forces force; + Vector6d previousVelocity{0}; Vector6d velocity{0}; double kineticEnergy{0}; Vector6d displacements{0};