manolas
/
MySources
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Added code for exporting pe.gradual application of external forces. average residual forces threshold

This commit is contained in:
iasonmanolas 2021-07-07 22:43:44 +03:00
parent 992f0fe6ae
commit 2509b0a795
6 changed files with 193 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

250
drmsimulationmodel.cpp
View File

@ -231,6 +231,7 @@ void DRMSimulationModel::reset()
Dt = mSettings.Dtini;
numOfDampings = 0;
shouldTemporarilyDampForces = false;
externalLoadStep = 1;
}
VectorType DRMSimulationModel::computeDisplacementDifferenceDerivative(
@ -958,7 +959,6 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
force.residual = force.external;
force.internal = 0;
}
double totalResidualForcesNorm = 0;
for (size_t ei = 0; ei < pMesh->EN(); ei++) {
for (int i = 0; i < 4; i++) {
std::pair<int, Vector6d> internalForcePair
@ -972,6 +972,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
force.residual = force.residual + (internalForcePair.second * -1);
}
}
double totalResidualForcesNorm = 0;
Vector6d sumOfResidualForces(0);
for (size_t vi = 0; vi < pMesh->VN(); vi++) {
Node::Forces &force = pMesh->nodes[vi].force;
@ -984,7 +985,8 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
}
pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm;
pMesh->totalResidualForcesNorm = totalResidualForcesNorm;
// mesh->totalResidualForcesNorm = sumOfResidualForces.norm();
pMesh->averageResidualForcesNorm = totalResidualForcesNorm / pMesh->VN();
// pMesh->averageResidualForcesNorm = sumOfResidualForces.norm() / pMesh->VN();
// plotYValues.push_back(totalResidualForcesNorm);
// auto xPlot = matplot::linspace(0, plotYValues.size(), plotYValues.size());
@ -1147,7 +1149,7 @@ DRMSimulationModel::DRMSimulationModel() {}
void DRMSimulationModel::updateNodalMasses()
{
double gamma = mSettings.gamma;
const size_t untilStep = 500;
const size_t untilStep = 4000;
if (shouldTemporarilyDampForces && mCurrentSimulationStep < untilStep) {
gamma *= 1e6 * (1 - static_cast<double>(mCurrentSimulationStep) / untilStep);
}
@ -1233,10 +1235,9 @@ void DRMSimulationModel::updateNodalVelocities()
const VertexIndex vi = pMesh->getIndex(v);
Node &node = pMesh->nodes[v];
node.velocity = node.velocity + node.acceleration * Dt;
// if (vi == 10) {
// std::cout << "Velocity:" << node.velocity[0] << " " << node.velocity[1] << " "
// << node.velocity[2] << std::endl;
// }
if (std::isnan(node.velocity.norm())) {
std::cout << "Velocity:" << node.velocity.toString() << std::endl;
}
}
updateKineticEnergy();
}
@ -1461,7 +1462,7 @@ void DRMSimulationModel::updateKineticEnergy()
+ node.mass.rotationalI3 * pow(node.velocity[4], 2)
+ node.mass.rotationalI2 * pow(node.velocity[5], 2));
node.kineticEnergy += nodeTranslationalKineticEnergy + nodeRotationalKineticEnergy;
node.kineticEnergy = nodeTranslationalKineticEnergy + nodeRotationalKineticEnergy;
assert(node.kineticEnergy < 1e15);
pMesh->currentTotalKineticEnergy += node.kineticEnergy;
@ -1959,26 +1960,10 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
const SimulationResults &solutionGuess)
{
assert(pJob->pMesh.operator bool());
auto t1 = std::chrono::high_resolution_clock::now();
auto beginTime = std::chrono::high_resolution_clock::now();
mSettings = settings;
reset();
double totalExternalForcesNorm = 0;
for (const std::pair<VertexIndex, Vector6d> &externalForce : pJob->nodalExternalForces) {
totalExternalForcesNorm += externalForce.second.norm();
}
if (!pJob->nodalExternalForces.empty()) {
mSettings.totalResidualForcesNormThreshold
= settings.totalExternalForcesNormPercentageTermination * totalExternalForcesNorm;
} 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;
}
history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
// if (!pJob->nodalExternalForces.empty()) {
@ -2044,12 +2029,38 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
}
updateNodalMasses();
updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
const size_t movingAverageSampleSize = 200;
std::queue<double> residualForcesMovingAverageHistorySample;
std::unordered_map<VertexIndex, Vector6d> 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(nodalExternalForces, constrainedVertices);
double averageExternalForcesNorm =
// sumOfExternalForces.norm() / pMesh->VN();
totalExternalForcesNorm / pMesh->VN();
if (!nodalExternalForces.empty()) {
mSettings.totalResidualForcesNormThreshold
= std::min(settings.totalExternalForcesNormPercentageTermination
* totalExternalForcesNorm,
1e-3);
} 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;
}
// const size_t movingAverageSampleSize = 200;
// std::queue<double> residualForcesMovingAverageHistorySample;
std::vector<double> percentageOfConvergence;
double residualForcesMovingAverageDerivativeMax = 0;
// double residualForcesMovingAverageDerivativeMax = 0;
while (settings.maxDRMIterations == 0 || mCurrentSimulationStep < settings.maxDRMIterations) {
if ((mSettings.isDebugMode && mCurrentSimulationStep == 50000)) {
// std::filesystem::create_directory("./PatternOptimizationNonConv");
@ -2087,38 +2098,46 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// }
updateElementalLengths();
mCurrentSimulationStep++;
double sumOfDisplacementsNorm = 0;
for (size_t vi = 0; vi < pMesh->VN(); vi++) {
sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm();
if (std::isnan(pMesh->currentTotalKineticEnergy)) {
std::cout << pMesh->currentTotalKineticEnergy << std::endl;
if (mSettings.beVerbose) {
std::cerr << "Infinite kinetic energy at step " << mCurrentSimulationStep
<< ". Exiting.." << std::endl;
}
break;
}
sumOfDisplacementsNorm /= pMesh->bbox.Diag();
pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm;
// pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm
// + mCurrentSimulationStep
// * pMesh->residualForcesMovingAverage)
// / (1 + mCurrentSimulationStep);
pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm
/ movingAverageSampleSize;
residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage);
if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) {
const double firstElementValue = residualForcesMovingAverageHistorySample.front();
pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize;
residualForcesMovingAverageHistorySample.pop();
pMesh->residualForcesMovingAverageDerivativeNorm
= std::abs((residualForcesMovingAverageHistorySample.back()
- residualForcesMovingAverageHistorySample.front()))
/ (movingAverageSampleSize);
residualForcesMovingAverageDerivativeMax
= std::max(pMesh->residualForcesMovingAverageDerivativeNorm,
residualForcesMovingAverageDerivativeMax);
pMesh->residualForcesMovingAverageDerivativeNorm
/= residualForcesMovingAverageDerivativeMax;
// std::cout << "Normalized derivative:"
// << residualForcesMovingAverageDerivativeNorm
// << std::endl;
}
// double sumOfDisplacementsNorm = 0;
// for (size_t vi = 0; vi < pMesh->VN(); vi++) {
// sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm();
// }
// sumOfDisplacementsNorm /= pMesh->bbox.Diag();
// pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm;
// // pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm
// // + mCurrentSimulationStep
// // * pMesh->residualForcesMovingAverage)
// // / (1 + mCurrentSimulationStep);
// pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm
// / movingAverageSampleSize;
// residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage);
// if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) {
// const double firstElementValue = residualForcesMovingAverageHistorySample.front();
// pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize;
// residualForcesMovingAverageHistorySample.pop();
// pMesh->residualForcesMovingAverageDerivativeNorm
// = std::abs((residualForcesMovingAverageHistorySample.back()
// - residualForcesMovingAverageHistorySample.front()))
// / (movingAverageSampleSize);
// residualForcesMovingAverageDerivativeMax
// = std::max(pMesh->residualForcesMovingAverageDerivativeNorm,
// residualForcesMovingAverageDerivativeMax);
// pMesh->residualForcesMovingAverageDerivativeNorm
// /= residualForcesMovingAverageDerivativeMax;
// // std::cout << "Normalized derivative:"
// // << residualForcesMovingAverageDerivativeNorm
// // << std::endl;
// }
// pMesh->previousTotalPotentialEnergykN =
// pMesh->currentTotalPotentialEnergykN;
@ -2127,10 +2146,21 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
if (mSettings.beVerbose && mSettings.isDebugMode
&& mCurrentSimulationStep % mSettings.debugModeStep == 0) {
printCurrentState();
std::cout << "Percentage of target:"
<< 100 * mSettings.totalResidualForcesNormThreshold
/ pMesh->totalResidualForcesNorm
<< "%" << std::endl;
// auto t2 = std::chrono::high_resolution_clock::now();
// const double executionTime_min
// = std::chrono::duration_cast<std::chrono::minutes>(t2 - beginTime).count();
// std::cout << "Execution time(min):" << executionTime_min << std::endl;
if (mSettings.useAverageResidualForcesNorm) {
const double percOfConv = averageExternalForcesNorm
/ pMesh->averageResidualForcesNorm;
std::cout << "Percentage of target:" << percOfConv << "%" << std::endl;
} else {
std::cout << "Percentage of target:"
<< 100 * mSettings.totalResidualForcesNormThreshold
/ pMesh->totalResidualForcesNorm
<< "%" << std::endl;
}
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver",
// movingAverages);
@ -2149,13 +2179,6 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
/ pMesh->totalResidualForcesNorm);
}
if (std::isnan(pMesh->currentTotalKineticEnergy)) {
if (mSettings.beVerbose) {
std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl;
}
break;
}
if (mSettings.shouldCreatePlots && mSettings.isDebugMode
&& mCurrentSimulationStep % mSettings.debugModeStep == 0) {
// SimulationResultsReporter::createPlot("Number of Steps",
@ -2164,18 +2187,18 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver",
// history.residualForcesMovingAverage);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Residual Forces",
// history.logResidualForces,
// {},
// history.redMarks);
SimulationResultsReporter::createPlot("Number of Steps",
"Log of Residual Forces",
history.logResidualForces,
{},
history.redMarks);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Kinetic energy",
// history.kineticEnergy,
// {},
// history.redMarks);
SimulationResultsReporter reporter;
reporter.reportHistory(history, "IntermediateResults", pJob->pMesh->getLabel());
// SimulationResultsReporter reporter;
// reporter.reportHistory(history, "IntermediateResults", pJob->pMesh->getLabel());
// SimulationResultsReporter::createPlot("Number of Iterations",
// "Sum of normalized displacement norms",
// history.sumOfNormalizedDisplacementNorms /*,
@ -2234,8 +2257,12 @@ mesh->currentTotalPotentialEnergykN*/
&& pMesh->currentTotalTranslationalKineticEnergy
< mSettings.totalTranslationalKineticEnergyThreshold
&& mCurrentSimulationStep > 20 && numOfDampings > 0;
const bool fullfillsAverageResidualForcesNormCriterion
= mSettings.useAverageResidualForcesNorm
&& pMesh->averageResidualForcesNorm / averageExternalForcesNorm < 1e-2;
const bool fullfillsResidualForcesNormTerminationCriterion
= pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold;
= mSettings.useResidualForcesNorm
&& pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold;
const bool fullfillsMovingAverageNormTerminationCriterion
= pMesh->residualForcesMovingAverage
< mSettings.residualForcesMovingAverageNormThreshold;
@ -2246,7 +2273,8 @@ mesh->currentTotalPotentialEnergykN*/
= fullfillsKineticEnergyTerminationCriterion
// || fullfillsMovingAverageNormTerminationCriterion
// || fullfillsMovingAverageDerivativeNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion;
|| fullfillsResidualForcesNormTerminationCriterion
|| fullfillsAverageResidualForcesNormCriterion;
if (shouldTerminate) {
if (mSettings.beVerbose && !mSettings.isDebugMode) {
std::cout << "Simulation converged." << std::endl;
@ -2265,29 +2293,52 @@ mesh->currentTotalPotentialEnergykN*/
<< std::endl;
}
}
break;
if (mSettings.desiredGradualExternalLoadsSteps == externalLoadStep) {
break;
} else {
externalLoadStep++;
std::unordered_map<VertexIndex, Vector6d> nodalExternalForces
= pJob->nodalExternalForces;
double totalExternalForcesNorm = 0;
for (auto &nodalForce : nodalExternalForces) {
const double percentageOfExternalLoads
= double(externalLoadStep) / mSettings.desiredGradualExternalLoadsSteps;
nodalForce.second = nodalForce.second * percentageOfExternalLoads;
totalExternalForcesNorm += nodalForce.second.norm();
}
updateNodalExternalForces(nodalExternalForces, constrainedVertices);
averageExternalForcesNorm = totalExternalForcesNorm / pMesh->VN();
if (!nodalExternalForces.empty()) {
mSettings.totalResidualForcesNormThreshold = 1e-2 * totalExternalForcesNorm;
}
if (mSettings.beVerbose) {
std::cout << "totalResidualForcesNormThreshold:"
<< mSettings.totalResidualForcesNormThreshold << std::endl;
}
}
// }
}
// const bool shouldCapDisplacements = mSettings.displacementCap.has_value();
// for (VertexType &v : pMesh->vert) {
// Node &node = pMesh->nodes[v];
// Vector6d stepDisplacement = node.velocity * Dt;
// if (shouldCapDisplacements
// && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) {
// stepDisplacement = stepDisplacement
// * (*mSettings.displacementCap
// / stepDisplacement.getTranslation().norm());
// }
// node.displacements = node.displacements - stepDisplacement;
// }
// applyDisplacements(constrainedVertices);
// if (!pJob->nodalForcedDisplacements.empty()) {
// applyForcedDisplacements(
const bool shouldCapDisplacements = mSettings.displacementCap.has_value();
for (VertexType &v : pMesh->vert) {
Node &node = pMesh->nodes[v];
Vector6d stepDisplacement = node.velocity * Dt;
if (shouldCapDisplacements
&& stepDisplacement.getTranslation().norm() > mSettings.displacementCap) {
stepDisplacement = stepDisplacement
* (*mSettings.displacementCap
/ stepDisplacement.getTranslation().norm());
}
node.displacements = node.displacements - stepDisplacement;
}
applyDisplacements(constrainedVertices);
if (!pJob->nodalForcedDisplacements.empty()) {
applyForcedDisplacements(
// pJob->nodalForcedDisplacements);
// }
// updateElementalLengths();
pJob->nodalForcedDisplacements);
}
updateElementalLengths();
// const double triggerPercentage = 0.01;
// const double xi_min = 0.55;
@ -2313,8 +2364,11 @@ mesh->currentTotalPotentialEnergykN*/
// draw();
}
}
auto endTime = std::chrono::high_resolution_clock::now();
SimulationResults results = computeResults(pJob);
results.executionTime = std::chrono::duration_cast<std::chrono::seconds>(endTime - beginTime)
.count();
if (mCurrentSimulationStep == settings.maxDRMIterations && mCurrentSimulationStep != 0) {
if (mSettings.beVerbose) {

4
drmsimulationmodel.hpp
View File

@ -39,10 +39,13 @@ public:
int maxDRMIterations{0};
bool shouldUseTranslationalKineticEnergyThreshold{false};
int gradualForcedDisplacementSteps{100};
int desiredGradualExternalLoadsSteps{1};
double gamma{0.8};
std::optional<double> displacementCap;
double totalResidualForcesNormThreshold{1e-3};
double totalExternalForcesNormPercentageTermination{1e-3};
bool useResidualForcesNorm{true};
bool useAverageResidualForcesNorm{false};
Settings() {}
};
@ -57,6 +60,7 @@ private:
matplot::line_handle plotHandle;
std::vector<double> plotYValues;
size_t numOfDampings{0};
int externalLoadStep{1};
const std::string meshPolyscopeLabel{"Simulation mesh"};
std::unique_ptr<SimulationMesh> pMesh;

24
reducedmodeloptimizer_structs.hpp
View File

@ -267,6 +267,7 @@ struct Colors
std::vector<double> perSimulationScenario_rotational;
std::vector<double> perSimulationScenario_total;
} objectiveValue;
std::vector<double> perScenario_fullPatternPotentialEnergy;
// std::vector<std::pair<std::string,double>> optimalXNameValuePairs;
std::vector<std::pair<std::string, double>> optimalXNameValuePairs;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
@ -284,6 +285,7 @@ struct Colors
inline static std::string Label{"Label"};
inline static std::string FullPatternLabel{"FullPatternLabel"};
inline static std::string Settings{"OptimizationSettings"};
inline static std::string FullPatternPotentialEnergies{"PE"};
};
void save(const string &saveToPath, const bool shouldExportDebugFiles = false)
@ -323,6 +325,16 @@ struct Colors
baseTriangle.cP2(0)[2]};
baseTriangleFullPattern.save(std::filesystem::path(saveToPath).string());
json_optimizationResults[JsonKeys::FullPatternLabel] = baseTriangleFullPattern.getLabel();
//potential energies
const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
std::vector<double> fullPatternPE(numberOfSimulationJobs);
for (int simulationScenarioIndex = 0; simulationScenarioIndex < numberOfSimulationJobs;
simulationScenarioIndex++) {
fullPatternPE[simulationScenarioIndex]
= perScenario_fullPatternPotentialEnergy[simulationScenarioIndex];
}
json_optimizationResults[JsonKeys::FullPatternPotentialEnergies] = fullPatternPE;
////Save to json file
std::filesystem::path jsonFilePath(
std::filesystem::path(saveToPath).append(JsonKeys::filename));
@ -652,6 +664,13 @@ struct Colors
os << "Obj value Rot " + simulationScenarioName;
os << "Obj value Total " + simulationScenarioName;
}
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < fullPatternSimulationJobs.size();
simulationScenarioIndex++) {
const std::string simulationScenarioName
= fullPatternSimulationJobs[simulationScenarioIndex]->getLabel();
os << "PE " + simulationScenarioName;
}
for (const auto &nameValuePair : optimalXNameValuePairs) {
os << nameValuePair.first;
}
@ -671,6 +690,11 @@ struct Colors
os << objectiveValue.perSimulationScenario_rotational[simulationScenarioIndex];
os << objectiveValue.perSimulationScenario_total[simulationScenarioIndex];
}
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < fullPatternSimulationJobs.size();
simulationScenarioIndex++) {
os << perScenario_fullPatternPotentialEnergy[simulationScenarioIndex];
}
for (const auto &optimalXNameValuePair : optimalXNameValuePairs) {
os << optimalXNameValuePair.second;
}

4
simulation_structs.hpp
View File

@ -492,7 +492,7 @@ struct SimulationResults
return labelPrefix + deliminator + job->pMesh->getLabel() + deliminator + job->getLabel();
}
void saveDeformedModel(const std::string &outputFolder = std::string())
bool saveDeformedModel(const std::string &outputFolder = std::string())
{
VCGEdgeMesh m;
vcg::tri::Append<VCGEdgeMesh, SimulationMesh>::MeshCopy(m, *job->pMesh);
@ -503,7 +503,7 @@ struct SimulationResults
displacements[vi][1],
displacements[vi][2]);
}
m.save(std::filesystem::path(outputFolder).append(getLabel() + ".ply").string());
return m.save(std::filesystem::path(outputFolder).append(getLabel() + ".ply").string());
}
void save(const std::string &outputFolder = std::string())

1
simulationmesh.hpp
View File

@ -43,6 +43,7 @@ public:
double currentTotalKineticEnergy{0};
double currentTotalTranslationalKineticEnergy{0};
double totalResidualForcesNorm{0};
double averageResidualForcesNorm{0};
double currentTotalPotentialEnergykN{0};
double previousTotalPotentialEnergykN{0};
double residualForcesMovingAverageDerivativeNorm{0};

10
utilities.hpp
View File

@ -116,6 +116,16 @@ struct Vector6d : public std::array<double, 6> {
{
return Eigen::Vector3d(this->operator[](3), this->operator[](4), this->operator[](5));
}
std::string toString() const
{
std::string s;
for (int i = 0; i < 6; i++) {
s.append(std::to_string(this->operator[](i)) + ",");
}
s.pop_back();
return s;
}
};
namespace Utilities {