MySources/drmsimulationmodel.hpp

301 lines
12 KiB
C++
Raw Normal View History

2020-11-27 11:47:21 +01:00
#ifndef BEAMFORMFINDER_HPP
#define BEAMFORMFINDER_HPP
//#ifdef USE_MATPLOT
2021-12-19 19:15:36 +01:00
//#include "matplot.h"
#include <matplot/matplot.h>
//#endif
2021-11-15 10:08:39 +01:00
#include "simulationmesh.hpp"
#include "simulationmodel.hpp"
2020-11-27 11:47:21 +01:00
#include <Eigen/Dense>
#include <filesystem>
#include <unordered_set>
#ifdef USE_ENSMALLEN
#include <armadillo>
#include <ensmallen.hpp>
#endif
2022-05-06 15:16:43 +02:00
class SimulationJob;
2020-11-27 11:47:21 +01:00
enum DoF { Ux = 0, Uy, Uz, Nx, Ny, Nr, NumDoF };
using DoFType = int;
using EdgeType = VCGEdgeMesh::EdgeType;
using VertexType = VCGEdgeMesh::VertexType;
struct DifferentiateWithRespectTo {
const VertexType &v;
const DoFType &dofi;
};
2021-01-12 13:40:25 +01:00
2021-11-15 10:08:39 +01:00
class DRMSimulationModel : public SimulationModel
2021-04-08 20:03:23 +02:00
{
2021-01-12 13:40:25 +01:00
public:
2021-05-24 13:43:32 +02:00
struct Settings
{
2022-05-06 15:16:43 +02:00
bool useTranslationalKineticEnergyForKineticDamping{true};
bool useTotalRotationalKineticEnergyForKineticDamping{false};
2021-05-24 13:43:32 +02:00
bool shouldDraw{false};
bool beVerbose{false};
bool shouldCreatePlots{false};
// int drawingStep{0};
// double residualForcesMovingAverageDerivativeNormThreshold{1e-8};
// double residualForcesMovingAverageNormThreshold{1e-8};
2021-05-24 13:43:32 +02:00
double Dtini{0.1};
double xi{0.9969};
2022-01-28 19:01:44 +01:00
std::optional<double> translationalKineticEnergyThreshold;
2021-07-28 17:38:05 +02:00
int gradualForcedDisplacementSteps{50};
// int desiredGradualExternalLoadsSteps{1};
double gamma{0.8};
2022-05-06 15:16:43 +02:00
double totalResidualForcesNormThreshold{1e-20};
2022-01-14 14:02:27 +01:00
double totalExternalForcesNormPercentageTermination{1e-5};
2021-07-28 17:38:05 +02:00
std::optional<int> maxDRMIterations;
std::optional<int> debugModeStep;
std::optional<double> totalTranslationalKineticEnergyThreshold;
std::optional<double> averageResidualForcesCriterionThreshold;
std::optional<double> linearGuessForceScaleFactor;
// std::optional<int> intermediateResultsSaveStep;
std::optional<bool> saveIntermediateBestStates;
2021-07-28 17:38:05 +02:00
std::optional<double> viscousDampingFactor;
2022-01-14 14:02:27 +01:00
Settings() {}
void save(const std::filesystem::path &jsonFilePath) const;
bool load(const std::filesystem::path &filePath);
struct JsonLabels
2021-07-28 17:38:05 +02:00
{
const std::string shouldDraw{"shouldDraw"};
const std::string beVerbose{"beVerbose"};
const std::string shouldCreatePlots{"shouldCreatePlots"};
const std::string Dtini{"DtIni"};
const std::string xi{"xi"};
2022-01-14 14:02:27 +01:00
const std::string gamma{"gamma"};
const std::string totalResidualForcesNormThreshold;
2021-07-28 17:38:05 +02:00
const std::string maxDRMIterations{"maxIterations"};
const std::string debugModeStep{"debugModeStep"};
const std::string totalTranslationalKineticEnergyThreshold{
"totalTranslationaKineticEnergyThreshold"};
const std::string averageResidualForcesCriterionThreshold{
"averageResidualForcesThreshold"};
const std::string linearGuessForceScaleFactor{"linearGuessForceScaleFactor"};
2022-01-14 14:02:27 +01:00
const std::string viscousDampingFactor{"viscousDampingFactor"};
2021-07-28 17:38:05 +02:00
};
static JsonLabels jsonLabels;
2021-05-24 13:43:32 +02:00
};
2020-11-27 11:47:21 +01:00
private:
Settings mSettings;
double Dt{mSettings.Dtini};
bool checkedForMaximumMoment{false};
bool shouldTemporarilyDampForces{false};
bool shouldTemporarilyAmplifyForces{true};
double externalMomentsNorm{0};
size_t mCurrentSimulationStep{0};
matplot::line_handle plotHandle;
std::vector<double> plotYValues;
size_t numOfDampings{0};
int externalLoadStep{1};
2021-07-15 11:01:12 +02:00
std::vector<bool> isVertexConstrained;
std::vector<bool> isRigidSupport;
double minTotalResidualForcesNorm{std::numeric_limits<double>::max()};
2020-11-27 11:47:21 +01:00
const std::string meshPolyscopeLabel{"Simulation mesh"};
std::unique_ptr<SimulationMesh> pMesh;
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> constrainedVertices;
SimulationHistory history;
// Eigen::Tensor<double, 4> theta3Derivatives;
// std::unordered_map<MyKeyType, double, key_hash> theta3Derivatives;
bool shouldApplyInitialDistortion = false;
//#ifdef USE_ENSMALLEN
// std::shared_ptr<SimulationJob> pJob;
//#endif
2021-11-15 10:08:39 +01:00
void reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings);
void updateNodalInternalForces(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void updateNodalExternalForces(
const std::unordered_map<VertexIndex, Vector6d> &nodalForces,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void updateResidualForces();
void updateRotationalDisplacements();
void updateElementalLengths();
2020-11-27 11:47:21 +01:00
void updateNodalMasses();
2020-11-27 11:47:21 +01:00
void updateNodalAccelerations();
2020-11-27 11:47:21 +01:00
void updateNodalVelocities();
2020-11-27 11:47:21 +01:00
void updateNodalDisplacements();
2020-11-27 11:47:21 +01:00
void applyForcedDisplacements(
const std::unordered_map<VertexIndex, Eigen::Vector3d> &nodalForcedDisplacements);
2020-11-27 11:47:21 +01:00
Vector6d computeElementTorsionalForce(const Element &element,
const Vector6d &displacementDifference,
const std::unordered_set<DoFType> &constrainedDof);
2020-11-27 11:47:21 +01:00
// BeamFormFinder::Vector6d computeElementFirstBendingForce(
// const Element &element, const Vector6d &displacementDifference,
// const std::unordered_set<gsl::index> &constrainedDof);
2020-11-27 11:47:21 +01:00
// BeamFormFinder::Vector6d computeElementSecondBendingForce(
// const Element &element, const Vector6d &displacementDifference,
// const std::unordered_set<gsl::index> &constrainedDof);
2020-11-27 11:47:21 +01:00
void updateKineticEnergy();
2020-11-27 11:47:21 +01:00
void resetVelocities();
SimulationResults computeResults(const std::shared_ptr<SimulationJob> &pJob);
void updateNodePosition(
VertexType &v,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void applyDisplacements(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
2020-11-27 11:47:21 +01:00
2021-04-30 12:13:58 +02:00
#ifdef POLYSCOPE_DEFINED
2022-05-06 15:16:43 +02:00
void draw(const std::shared_ptr<SimulationJob> &pJob, const std::string &screenshotsFolder = {});
2021-04-30 12:13:58 +02:00
#endif
2020-11-27 11:47:21 +01:00
void
updateNodalInternalForce(Vector6d &nodalInternalForce,
const Vector6d &elementInternalForce,
const std::unordered_set<DoFType> &nodalFixedDof);
2020-12-03 19:56:03 +01:00
Vector6d computeElementInternalForce(
2020-11-27 11:47:21 +01:00
const Element &elem, const Node &n0, const Node &n1,
const std::unordered_set<DoFType> &n0ConstrainedDof,
const std::unordered_set<DoFType> &n1ConstrainedDof);
2020-12-03 19:56:03 +01:00
Vector6d computeElementAxialForce(const ::EdgeType &e) const;
2020-11-27 11:47:21 +01:00
VectorType computeDisplacementDifferenceDerivative(
const EdgeType &e, const DifferentiateWithRespectTo &dui) const;
double
computeDerivativeElementLength(const EdgeType &e,
const DifferentiateWithRespectTo &dui) const;
VectorType computeDerivativeT1(const EdgeType &e,
const DifferentiateWithRespectTo &dui) const;
VectorType
computeDerivativeOfNormal(const VertexType &v,
const DifferentiateWithRespectTo &dui) const;
VectorType computeDerivativeT3(const EdgeType &e,
const DifferentiateWithRespectTo &dui) const;
VectorType computeDerivativeT2(const EdgeType &e,
const DifferentiateWithRespectTo &dui) const;
double computeDerivativeTheta2(const EdgeType &e, const VertexIndex &evi,
const VertexIndex &dwrt_evi,
const DoFType &dwrt_dofi) const;
void updateElementalFrames();
2021-07-21 16:59:46 +02:00
VectorType computeDerivativeOfR(const EdgeType &e,
const DifferentiateWithRespectTo &dui) const;
2020-11-27 11:47:21 +01:00
static double computeDerivativeOfNorm(const VectorType &x,
const VectorType &derivativeOfX);
static VectorType computeDerivativeOfCrossProduct(
const VectorType &a, const VectorType &derivativeOfA, const VectorType &b,
const VectorType &derivativeOfB);
double computeTheta3(const EdgeType &e, const VertexType &v);
double computeDerivativeTheta3(const EdgeType &e, const VertexType &v,
const DifferentiateWithRespectTo &dui) const;
2021-04-08 20:03:23 +02:00
double computeTotalPotentialEnergy();
2020-11-27 11:47:21 +01:00
void computeRigidSupports();
void updateNormalDerivatives();
void updateT1Derivatives();
void updateT2Derivatives();
void updateT3Derivatives();
void updateRDerivatives();
double computeDerivativeTheta1(const EdgeType &e, const VertexIndex &evi,
const VertexIndex &dwrt_evi,
const DoFType &dwrt_dofi) const;
// void updatePositionsOnTheFly(
// const std::unordered_map<VertexIndex,
// std::unordered_set<gsl::index>>
// &fixedVertices);
void updateResidualForcesOnTheFly(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices);
void updatePositionsOnTheFly(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices);
2020-12-03 19:56:03 +01:00
void updateNodeNormal(
VertexType &v,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices);
2020-12-09 16:59:18 +01:00
void applyForcedNormals(
const std::unordered_map<VertexIndex, VectorType> nodalForcedRotations);
2021-01-14 14:39:19 +01:00
void printCurrentState() const;
2021-01-12 13:40:25 +01:00
void printDebugInfo() const;
2021-04-08 20:03:23 +02:00
double computeTotalInternalPotentialEnergy();
void applySolutionGuess(const SimulationResults &solutionGuess,
const std::shared_ptr<SimulationJob> &pJob);
void updateNodeNr(VertexType &v);
2021-07-28 17:38:05 +02:00
void applyKineticDamping(const std::shared_ptr<SimulationJob> &pJob);
2021-11-15 10:08:39 +01:00
virtual SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob);
void reset(const std::shared_ptr<SimulationJob> &pJob);
std::vector<std::array<Vector6d, 4>> computeInternalForces(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
2021-04-08 20:03:23 +02:00
public:
DRMSimulationModel();
SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
2021-11-15 10:08:39 +01:00
const Settings &settings,
const SimulationResults &solutionGuess = SimulationResults());
2022-05-06 15:16:43 +02:00
//#ifdef USE_ENSMALLEN
// std::shared_ptr<SimulationJob> pJob;
// double EvaluateWithGradient(const arma::mat &x, arma::mat &g);
// void setJob(const std::shared_ptr<SimulationJob> &pJob);
// SimulationMesh *getDeformedMesh(const arma::mat &x, const std::shared_ptr<SimulationJob> &pJob);
// double Evaluate(const arma::mat &x);
//#endif
2020-12-09 16:59:18 +01:00
2021-05-24 13:43:32 +02:00
static void runUnitTests();
2020-11-27 11:47:21 +01:00
};
2021-07-28 17:38:05 +02:00
inline DRMSimulationModel::Settings::JsonLabels DRMSimulationModel::Settings::jsonLabels;
2020-11-27 11:47:21 +01:00
template <typename PointType> PointType Cross(PointType p1, PointType p2) {
return p1 ^ p2;
}
inline size_t currentStep{0};
inline bool TriggerBreakpoint(const VertexIndex &vi, const EdgeIndex &ei,
const DoFType &dofi) {
const size_t numberOfVertices = 10;
const VertexIndex middleNodeIndex = numberOfVertices / 2;
// return vi == middleNodeIndex && dofi == 1;
return dofi == 1 && ((vi == 1 && ei == 0) || (vi == 9 && ei == 9));
}
inline bool TriggerBreakpoint(const VertexIndex &vi, const EdgeIndex &ei) {
const size_t numberOfVertices = 10;
const VertexIndex middleNodeIndex = numberOfVertices / 2;
return (vi == middleNodeIndex);
// return (vi == 0 || vi == numberOfVertices - 1) && currentStep == 1;
return (vi == 1 && ei == 0) || (vi == 9 && ei == 9);
}
#endif // BEAMFORMFINDER_HPP