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

Using quadratic reset. Unit tests 2120 ,3333,9268 iterations

This commit is contained in:
iasonmanolas 2022-07-26 20:04:07 +03:00
parent 76fb1b8609
commit 90551e5485
4 changed files with 209 additions and 185 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

238
drmsimulationmodel.cpp
View File

@ -45,10 +45,11 @@ void DRMSimulationModel::runUnitTests()
Settings settings; Settings settings;
settings.Dtini = 0.1; settings.Dtini = 0.1;
settings.xi = 0.9969; settings.xi = 0.9969;
settings.maxDRMIterations = 0; settings.totalResidualForcesNormThreshold = 1;
formFinder.mSettings.totalResidualForcesNormThreshold = 1;
settings.shouldDraw = false; settings.shouldDraw = false;
settings.beVerbose = true; settings.beVerbose = true;
// settings.debugModeStep = 1000;
// settings.shouldDraw = true;
settings.shouldCreatePlots = true; settings.shouldCreatePlots = true;
SimulationResults simpleBeam_simulationResults SimulationResults simpleBeam_simulationResults
= formFinder.executeSimulation(std::make_shared<SimulationJob>(beamSimulationJob), settings); = formFinder.executeSimulation(std::make_shared<SimulationJob>(beamSimulationJob), settings);
@ -104,10 +105,11 @@ void DRMSimulationModel::runUnitTests()
shortSpanGridshellSimulationJob.pMesh->setBeamMaterial(0.3, 200 * 1e9); shortSpanGridshellSimulationJob.pMesh->setBeamMaterial(0.3, 200 * 1e9);
assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
shortSpanGridshellSimulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026}); shortSpanGridshellSimulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026});
DRMSimulationModel formFinder2;
SimulationResults shortSpanGridshellSimulationResults SimulationResults shortSpanGridshellSimulationResults
= formFinder.executeSimulation(std::make_shared<SimulationJob>( = formFinder2.executeSimulation(std::make_shared<SimulationJob>(
shortSpanGridshellSimulationJob), shortSpanGridshellSimulationJob),
settings); settings);
shortSpanGridshellSimulationResults.save(); shortSpanGridshellSimulationResults.save();
const std::string groundOfTruthBinaryFilename const std::string groundOfTruthBinaryFilename
@ -183,10 +185,11 @@ void DRMSimulationModel::runUnitTests()
} }
longSpanGridshell_simulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026}); longSpanGridshell_simulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026});
DRMSimulationModel formFinder3;
SimulationResults longSpanGridshell_simulationResults SimulationResults longSpanGridshell_simulationResults
= formFinder.executeSimulation(std::make_shared<SimulationJob>( = formFinder3.executeSimulation(std::make_shared<SimulationJob>(
longSpanGridshell_simulationJob), longSpanGridshell_simulationJob),
settings); settings);
longSpanGridshell_simulationResults.save(); longSpanGridshell_simulationResults.save();
const std::string longSpan_groundOfTruthBinaryFilename const std::string longSpan_groundOfTruthBinaryFilename
@ -217,22 +220,14 @@ void DRMSimulationModel::runUnitTests()
// polyscope::show(); // polyscope::show();
} }
void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob) void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob)
{ {
mCurrentSimulationStep = 0; //#ifdef USE_ENSMALLEN
history.clear(); // this->pJob = pJob;
history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel(); //#endif
plotYValues.clear(); pMesh.reset();
plotHandle.reset(); pMesh = std::make_unique<SimulationMesh>(*pJob->pMesh);
checkedForMaximumMoment = false; vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
externalMomentsNorm = 0;
Dt = mSettings.Dtini;
numOfDampings = 0;
shouldTemporarilyDampForces = false;
externalLoadStep = 1;
isVertexConstrained.clear();
minTotalResidualForcesNorm = std::numeric_limits<double>::max();
constrainedVertices.clear(); constrainedVertices.clear();
constrainedVertices = pJob->constrainedVertices; constrainedVertices = pJob->constrainedVertices;
@ -247,6 +242,26 @@ void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob)
for (auto fixedVertex : constrainedVertices) { for (auto fixedVertex : constrainedVertices) {
isVertexConstrained[fixedVertex.first] = true; isVertexConstrained[fixedVertex.first] = true;
} }
}
void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &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<double>::max();
reset(pJob);
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) { if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) {
@ -257,14 +272,14 @@ void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob)
polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel + "_" + pMesh->getLabel(), polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel + "_" + pMesh->getLabel(),
pMesh->getEigenVertices(), pMesh->getEigenVertices(),
pMesh->getEigenEdges()) pMesh->getEigenEdges())
->setRadius(pMesh->elements[0].dimensions.getDrawingRadius()); ->setRadius(0.002);
// registerWorldAxes(); // registerWorldAxes();
} }
#endif #endif
// if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) { if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) {
// shouldApplyInitialDistortion = true; shouldApplyInitialDistortion = true;
// } }
updateElementalFrames(); updateElementalFrames();
} }
@ -1414,7 +1429,6 @@ void DRMSimulationModel::updateResidualForces()
void DRMSimulationModel::computeRigidSupports() void DRMSimulationModel::computeRigidSupports()
{ {
assert(pMesh != nullptr);
isRigidSupport.clear(); isRigidSupport.clear();
isRigidSupport.resize(pMesh->VN(), false); isRigidSupport.resize(pMesh->VN(), false);
for (const VertexType &v : pMesh->vert) { for (const VertexType &v : pMesh->vert) {
@ -1625,6 +1639,7 @@ void DRMSimulationModel::updateNodalVelocities()
for (VertexType &v : pMesh->vert) { for (VertexType &v : pMesh->vert) {
const VertexIndex vi = pMesh->getIndex(v); const VertexIndex vi = pMesh->getIndex(v);
Node &node = pMesh->nodes[v]; Node &node = pMesh->nodes[v];
node.previousVelocity = node.velocity;
if (mSettings.viscousDampingFactor.has_value()) { if (mSettings.viscousDampingFactor.has_value()) {
const Vector6d massOverDt = node.mass_6d / Dt; const Vector6d massOverDt = node.mass_6d / Dt;
// const Vector6d visciousDampingFactor(viscuousDampingConstant / 2); // const Vector6d visciousDampingFactor(viscuousDampingConstant / 2);
@ -1777,11 +1792,17 @@ void DRMSimulationModel::updateNodeNormal(
void DRMSimulationModel::applyDisplacements( void DRMSimulationModel::applyDisplacements(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices) const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices)
{ {
for (VertexType &v : pMesh->vert) { std::for_each(
updateNodePosition(v, fixedVertices); #ifdef ENABLE_PARALLEL_DRM
updateNodeNormal(v, fixedVertices); std::execution::par_unseq,
updateNodeNr(v); #endif
} pMesh->vert.begin(),
pMesh->vert.end(),
[&](auto &v) {
updateNodePosition(v, fixedVertices);
updateNodeNormal(v, fixedVertices);
updateNodeNr(v);
});
updateElementalFrames(); updateElementalFrames();
if (mSettings.shouldDraw) { if (mSettings.shouldDraw) {
pMesh->updateEigenEdgeAndVertices(); pMesh->updateEigenEdgeAndVertices();
@ -1847,8 +1868,12 @@ void DRMSimulationModel::applyForcedNormals(
// NOTE: Is this correct? Should the kinetic energy be computed like that? // NOTE: Is this correct? Should the kinetic energy be computed like that?
void DRMSimulationModel::updateKineticEnergy() void DRMSimulationModel::updateKineticEnergy()
{ {
pMesh->pre_previousTotalKineticEnergy = pMesh->previousTotalKineticEnergy;
pMesh->pre_previousTotalTranslationalKineticEnergy
= pMesh->previousTotalTranslationalKineticEnergy;
pMesh->pre_previousTotalRotationalKineticEnergy = pMesh->previousTotalRotationalKineticEnergy;
pMesh->previousTotalKineticEnergy = pMesh->currentTotalKineticEnergy; pMesh->previousTotalKineticEnergy = pMesh->currentTotalKineticEnergy;
pMesh->previousTranslationalKineticEnergy = pMesh->currentTotalTranslationalKineticEnergy; pMesh->previousTotalTranslationalKineticEnergy = pMesh->currentTotalTranslationalKineticEnergy;
pMesh->previousTotalRotationalKineticEnergy = pMesh->currentTotalRotationalKineticEnergy; pMesh->previousTotalRotationalKineticEnergy = pMesh->currentTotalRotationalKineticEnergy;
pMesh->currentTotalKineticEnergy = 0; pMesh->currentTotalKineticEnergy = 0;
pMesh->currentTotalTranslationalKineticEnergy = 0; pMesh->currentTotalTranslationalKineticEnergy = 0;
@ -1889,7 +1914,14 @@ void DRMSimulationModel::resetVelocities()
// // reset // // reset
// // current to 0 or this? // // current to 0 or this?
0; 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(); updateKineticEnergy();
} }
@ -1960,6 +1992,7 @@ void DRMSimulationModel::updatePositionsOnTheFly(
for (VertexType &v : pMesh->vert) { for (VertexType &v : pMesh->vert) {
Node &node = pMesh->nodes[v]; Node &node = pMesh->nodes[v];
node.previousVelocity = node.velocity;
node.velocity = node.velocity + node.acceleration * Dt; node.velocity = node.velocity + node.acceleration * Dt;
} }
updateKineticEnergy(); updateKineticEnergy();
@ -2564,15 +2597,26 @@ void DRMSimulationModel::applyKineticDamping(const std::shared_ptr<SimulationJob
{ {
// if (!mSettings.viscousDampingFactor.has_value()) { // if (!mSettings.viscousDampingFactor.has_value()) {
// const bool shouldCapDisplacements = mSettings.displacementCap.has_value(); // const bool shouldCapDisplacements = mSettings.displacementCap.has_value();
// const double &KE1 = pMesh->pre_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) { for (VertexType &v : pMesh->vert) {
Node &node = pMesh->nodes[v]; Node &node = pMesh->nodes[v];
Vector6d stepDisplacement = node.velocity * 0.5 * Dt; // Vector6d stepDisplacement = node.velocity * 0.5 * Dt;
// if (shouldCapDisplacements // if (shouldCapDisplacements
// && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) { // && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) {
// stepDisplacement = stepDisplacement // stepDisplacement = stepDisplacement
// * (*mSettings.displacementCap // * (*mSettings.displacementCap
// / stepDisplacement.getTranslation().norm()); // / stepDisplacement.getTranslation().norm());
// } // }
// Vector6d stepDisplacement = node.velocity * 0.5 * Dt;
Vector6d stepDisplacement = node.velocity * Dt + node.previousVelocity * bitaDt;
node.displacements = node.displacements - stepDisplacement; node.displacements = node.displacements - stepDisplacement;
} }
applyDisplacements(constrainedVertices); applyDisplacements(constrainedVertices);
@ -2598,51 +2642,38 @@ void DRMSimulationModel::applyKineticDamping(const std::shared_ptr<SimulationJob
++numOfDampings; ++numOfDampings;
} }
void DRMSimulationModel::updateDerivatives()
{
updateNormalDerivatives();
updateT1Derivatives();
updateRDerivatives();
updateT2Derivatives();
updateT3Derivatives();
}
SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob) SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob)
{ {
assert(pMesh != nullptr && pMesh->VN() == pJob->pMesh->VN() && pMesh->EN() == pJob->pMesh->EN());
reset(pJob);
auto beginTime = std::chrono::high_resolution_clock::now(); auto beginTime = std::chrono::high_resolution_clock::now();
updateNodalMasses(); 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<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(pJob->nodalExternalForces, constrainedVertices); updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
if (!pJob->nodalExternalForces.empty()) { if (!pJob->nodalExternalForces.empty()
&& !mSettings.totalResidualForcesNormThreshold.has_value()) {
mSettings.totalResidualForcesNormThreshold mSettings.totalResidualForcesNormThreshold
= mSettings.totalExternalForcesNormPercentageTermination = mSettings.totalExternalForcesNormPercentageTermination
* pMesh->totalExternalForcesNorm; * pMesh->totalExternalForcesNorm;
} else { } /*else {
mSettings.totalResidualForcesNormThreshold = 1e-3; mSettings.totalResidualForcesNormThreshold = 1e-3;
std::cout << "No forces setted default residual forces norm threshold" << std::endl; std::cout << "No forces setted default residual forces norm threshold" << std::endl;
} }*/
if (mSettings.beVerbose) { if (mSettings.beVerbose) {
std::cout << "totalResidualForcesNormThreshold:" // std::cout << "totalResidualForcesNormThreshold:"
<< mSettings.totalResidualForcesNormThreshold << std::endl; // << mSettings.totalResidualForcesNormThreshold << std::endl;
if (mSettings.averageResidualForcesCriterionThreshold.has_value()) { if (mSettings.averageResidualForcesCriterionThreshold.has_value()) {
std::cout << "average/extNorm threshold:" std::cout << "average/extNorm threshold:"
<< *mSettings.averageResidualForcesCriterionThreshold << std::endl; << *mSettings.averageResidualForcesCriterionThreshold << std::endl;
}
} }
}
if (mSettings.beVerbose) { if (mSettings.beVerbose) {
std::cout << "Executing simulation for mesh with:" << pMesh->VN() << " nodes and " std::cout << "Executing simulation for mesh with:" << pMesh->VN() << " nodes and "
@ -2655,25 +2686,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// double residualForcesMovingAverageDerivativeMax = 0; // double residualForcesMovingAverageDerivativeMax = 0;
while (!mSettings.maxDRMIterations.has_value() while (!mSettings.maxDRMIterations.has_value()
|| mCurrentSimulationStep < mSettings.maxDRMIterations.value()) { || mCurrentSimulationStep < mSettings.maxDRMIterations.value()) {
// if ((mSettings.debugModeStep.has_value() && mCurrentSimulationStep == 50000)) { updateDerivatives();
// std::filesystem::create_directory("./PatternOptimizationNonConv");
// pJob->save("./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++;
}
updateResidualForcesOnTheFly(constrainedVertices); updateResidualForcesOnTheFly(constrainedVertices);
// TODO: write parallel function for updating positions // TODO: write parallel function for updating positions
@ -2690,9 +2703,6 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
pJob->nodalForcedDisplacements); pJob->nodalForcedDisplacements);
} }
// if (!pJob->nodalForcedNormals.empty()) {
// applyForcedNormals(pJob->nodalForcedNormals);
// }
updateElementalLengths(); updateElementalLengths();
mCurrentSimulationStep++; mCurrentSimulationStep++;
if (std::isnan(pMesh->currentTotalKineticEnergy)) { if (std::isnan(pMesh->currentTotalKineticEnergy)) {
@ -2790,8 +2800,8 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// << std::endl; // << std::endl;
} }
history.stepPulse(*pMesh); history.stepPulse(*pMesh);
percentageOfConvergence.push_back(100 * mSettings.totalResidualForcesNormThreshold // percentageOfConvergence.push_back(100 * mSettings.totalResidualForcesNormThreshold
/ pMesh->totalResidualForcesNorm); // / pMesh->totalResidualForcesNorm);
} }
if (mSettings.shouldCreatePlots && mSettings.debugModeStep.has_value() if (mSettings.shouldCreatePlots && mSettings.debugModeStep.has_value()
@ -2880,15 +2890,19 @@ currentSimulationStep > maxDRMIterations*/
// std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm
// << std::endl; // << std::endl;
const bool fullfillsResidualForcesNormTerminationCriterion const bool fullfillsResidualForcesNormTerminationCriterion
= !mSettings.averageResidualForcesCriterionThreshold.has_value() = mSettings.averageResidualForcesCriterionThreshold.has_value()
&& pMesh->totalResidualForcesNorm / pMesh->totalExternalForcesNorm && pMesh->totalResidualForcesNorm / pMesh->totalExternalForcesNorm
< mSettings.totalExternalForcesNormPercentageTermination; < mSettings.totalExternalForcesNormPercentageTermination;
const bool fullfillsAverageResidualForcesNormTerminationCriterion const bool fullfillsAverageResidualForcesNormTerminationCriterion
= mSettings.averageResidualForcesCriterionThreshold.has_value() = mSettings.averageResidualForcesCriterionThreshold.has_value()
&& (pMesh->totalResidualForcesNorm / pMesh->VN()) / pMesh->totalExternalForcesNorm && (pMesh->totalResidualForcesNorm / pMesh->VN()) / pMesh->totalExternalForcesNorm
< mSettings.averageResidualForcesCriterionThreshold.value(); < mSettings.averageResidualForcesCriterionThreshold.value();
const bool fullfillsResidualForcesNormThreshold
= mSettings.totalResidualForcesNormThreshold.has_value()
&& pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold;
if ((fullfillsAverageResidualForcesNormTerminationCriterion if ((fullfillsAverageResidualForcesNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion) || fullfillsResidualForcesNormTerminationCriterion
|| fullfillsResidualForcesNormThreshold)
&& numOfDampings > 0 && numOfDampings > 0
&& (pJob->nodalForcedDisplacements.empty() && (pJob->nodalForcedDisplacements.empty()
|| mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) { || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) {
@ -2906,12 +2920,15 @@ currentSimulationStep > maxDRMIterations*/
break; break;
} }
const bool isKineticEnergyPeak = (mSettings.useTotalRotationalKineticEnergyForKineticDamping const bool isKineticEnergyPeak
&& pMesh->previousTotalRotationalKineticEnergy = /* pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy
> pMesh->currentTotalRotationalKineticEnergy) ||*/
|| (mSettings.useTranslationalKineticEnergyForKineticDamping (mSettings.useTotalRotationalKineticEnergyForKineticDamping
&& pMesh->previousTranslationalKineticEnergy && pMesh->previousTotalRotationalKineticEnergy
> pMesh->currentTotalTranslationalKineticEnergy); > pMesh->currentTotalRotationalKineticEnergy)
|| (mSettings.useTranslationalKineticEnergyForKineticDamping
&& pMesh->previousTotalTranslationalKineticEnergy
> pMesh->currentTotalTranslationalKineticEnergy);
if (isKineticEnergyPeak) { if (isKineticEnergyPeak) {
const bool fullfillsKineticEnergyTerminationCriterion const bool fullfillsKineticEnergyTerminationCriterion
= mSettings.totalTranslationalKineticEnergyThreshold.has_value() = mSettings.totalTranslationalKineticEnergyThreshold.has_value()
@ -2964,13 +2981,13 @@ currentSimulationStep > maxDRMIterations*/
// } // }
} }
if (mSettings.useTranslationalKineticEnergyForKineticDamping) { // if (mSettings.useTranslationalKineticEnergyForKineticDamping) {
applyKineticDamping(pJob); applyKineticDamping(pJob);
if (mSettings.shouldCreatePlots) {
history.redMarks.push_back(mCurrentSimulationStep);
}
}
Dt *= mSettings.xi; Dt *= mSettings.xi;
if (mSettings.shouldCreatePlots) {
history.redMarks.push_back(mCurrentSimulationStep);
}
// }
// if (mSettings.isDebugMode) { // if (mSettings.isDebugMode) {
// std::cout << Dt << std::endl; // std::cout << Dt << std::endl;
// } // }
@ -3002,11 +3019,9 @@ currentSimulationStep > maxDRMIterations*/
void DRMSimulationModel::setStructure(const std::shared_ptr<SimulationMesh> &pMesh) void DRMSimulationModel::setStructure(const std::shared_ptr<SimulationMesh> &pMesh)
{ {
this->pMesh.reset(); std::cout << "This function is currently not implemented" << std::endl;
this->pMesh = std::make_unique<SimulationMesh>(*pMesh); assert(false);
vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh); std::terminate();
// assert(false);
// std::terminate();
} }
SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob, SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
@ -3014,8 +3029,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
const SimulationResults &solutionGuess) const SimulationResults &solutionGuess)
{ {
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now(); std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
mSettings = settings; reset(pJob, settings);
setStructure(pJob->pMesh);
assert(pJob->pMesh != nullptr); assert(pJob->pMesh != nullptr);
if (!solutionGuess.displacements.empty()) { if (!solutionGuess.displacements.empty()) {

9
drmsimulationmodel.hpp
View File

@ -47,8 +47,8 @@ public:
int gradualForcedDisplacementSteps{50}; int gradualForcedDisplacementSteps{50};
// int desiredGradualExternalLoadsSteps{1}; // int desiredGradualExternalLoadsSteps{1};
double gamma{0.8}; double gamma{0.8};
double totalResidualForcesNormThreshold{1e-20}; std::optional<double> totalResidualForcesNormThreshold;
double totalExternalForcesNormPercentageTermination{1e-3}; double totalExternalForcesNormPercentageTermination{1e-5};
std::optional<int> maxDRMIterations; std::optional<int> maxDRMIterations;
std::optional<int> debugModeStep; std::optional<int> debugModeStep;
std::optional<double> totalTranslationalKineticEnergyThreshold; std::optional<double> totalTranslationalKineticEnergyThreshold;
@ -104,10 +104,11 @@ private:
SimulationHistory history; SimulationHistory history;
// Eigen::Tensor<double, 4> theta3Derivatives; // Eigen::Tensor<double, 4> theta3Derivatives;
// std::unordered_map<MyKeyType, double, key_hash> theta3Derivatives; // std::unordered_map<MyKeyType, double, key_hash> theta3Derivatives;
bool shouldApplyInitialDistortion = false; bool shouldApplyInitialDistortion{false};
//#ifdef USE_ENSMALLEN //#ifdef USE_ENSMALLEN
// std::shared_ptr<SimulationJob> pJob; // std::shared_ptr<SimulationJob> pJob;
//#endif //#endif
void reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings);
void updateNodalInternalForces( void updateNodalInternalForces(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices); const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void updateNodalExternalForces( void updateNodalExternalForces(
@ -259,6 +260,8 @@ private:
std::vector<std::array<Vector6d, 4>> computeInternalForces( std::vector<std::array<Vector6d, 4>> computeInternalForces(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices); const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void updateDerivatives();
public: public:
DRMSimulationModel(); DRMSimulationModel();
SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob, SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob,

12
simulationmesh.cpp
View File

@ -268,12 +268,12 @@ void SimulationMesh::unregister() const
} }
#endif #endif
void SimulationMesh::setBeamCrossSection( void SimulationMesh::setBeamCrossSection(const CrossSectionType &beamDimensions)
const CrossSectionType &beamDimensions) { {
for (size_t ei = 0; ei < EN(); ei++) { for (size_t ei = 0; ei < EN(); ei++) {
elements[ei].dimensions = beamDimensions; elements[ei].dimensions = beamDimensions;
elements[ei].updateRigidity(); elements[ei].updateRigidity();
} }
} }
void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) { void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) {

135
simulationmesh.hpp
View File

@ -8,8 +8,8 @@
//extern bool drawGlobal; //extern bool drawGlobal;
struct Element; struct Element;
struct Node; struct Node;
using CrossSectionType = RectangularBeamDimensions; //using CrossSectionType = RectangularBeamDimensions;
//using CrossSectionType = CylindricalBeamDimensions; using CrossSectionType = CylindricalBeamDimensions;
class SimulationMesh : public VCGEdgeMesh { class SimulationMesh : public VCGEdgeMesh {
private: private:
@ -38,8 +38,11 @@ public:
std::vector<CrossSectionType> getBeamDimensions(); std::vector<CrossSectionType> getBeamDimensions();
std::vector<ElementMaterial> getBeamMaterial(); std::vector<ElementMaterial> getBeamMaterial();
double pre_previousTotalKineticEnergy{0};
double pre_previousTotalTranslationalKineticEnergy{0};
double pre_previousTotalRotationalKineticEnergy{0};
double previousTotalKineticEnergy{0}; double previousTotalKineticEnergy{0};
double previousTranslationalKineticEnergy{0}; double previousTotalTranslationalKineticEnergy{0};
double previousTotalRotationalKineticEnergy{0}; double previousTotalRotationalKineticEnergy{0};
double previousTotalResidualForcesNorm{0}; double previousTotalResidualForcesNorm{0};
double currentTotalKineticEnergy{0}; double currentTotalKineticEnergy{0};
@ -68,73 +71,76 @@ public:
#endif #endif
}; };
struct Element { struct Element
{
CrossSectionType dimensions;
ElementMaterial material;
CrossSectionType dimensions; void computeMaterialProperties(const ElementMaterial &material);
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); struct LocalFrame
// void computeDimensionsProperties(const RectangularBeamDimensions &dimensions); {
// void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions); VectorType t1;
void setDimensions(const CrossSectionType &dimensions); VectorType t2;
void setMaterial(const ElementMaterial &material); VectorType t3;
double getMass(const double &matrialDensity); };
struct LocalFrame { EdgeIndex ei;
VectorType t1; double length{0};
VectorType t2; double initialLength;
VectorType t3; LocalFrame frame;
};
EdgeIndex ei; struct Rigidity
double length{0}; {
double initialLength; double axial;
LocalFrame frame; 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 { VectorType f1_j;
double axial; VectorType f1_jplus1;
double torsional; VectorType f2_j;
double firstBending; VectorType f2_jplus1;
double secondBending; VectorType f3_j;
std::string toString() const { VectorType f3_jplus1;
return std::string("Rigidity:") + std::string("\nAxial=") + double cosRotationAngle_j;
std::to_string(axial) + std::string("\nTorsional=") + double cosRotationAngle_jplus1;
std::to_string(torsional) + std::string("\nFirstBending=") + double sinRotationAngle_j;
std::to_string(firstBending) + std::string("\nSecondBending=") + double sinRotationAngle_jplus1;
std::to_string(secondBending); std::vector<std::vector<VectorType>> derivativeT1;
} std::vector<std::vector<VectorType>> derivativeT2;
}; std::vector<std::vector<VectorType>> derivativeT3;
Rigidity rigidity; std::vector<VectorType> derivativeT1_j;
void updateRigidity(); std::vector<VectorType> derivativeT1_jplus1;
std::vector<VectorType> derivativeT2_j;
std::vector<VectorType> derivativeT2_jplus1;
std::vector<VectorType> derivativeT3_j;
std::vector<VectorType> derivativeT3_jplus1;
std::vector<VectorType> derivativeR_j;
std::vector<VectorType> derivativeR_jplus1;
struct RotationalDisplacements
{
double theta1{0}, theta2{0}, theta3{0};
};
RotationalDisplacements rotationalDisplacements_j;
RotationalDisplacements rotationalDisplacements_jplus1;
VectorType f1_j; static void computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A);
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<std::vector<VectorType>> derivativeT1;
std::vector<std::vector<VectorType>> derivativeT2;
std::vector<std::vector<VectorType>> derivativeT3;
std::vector<VectorType> derivativeT1_j;
std::vector<VectorType> derivativeT1_jplus1;
std::vector<VectorType> derivativeT2_j;
std::vector<VectorType> derivativeT2_jplus1;
std::vector<VectorType> derivativeT3_j;
std::vector<VectorType> derivativeT3_jplus1;
std::vector<VectorType> derivativeR_j;
std::vector<VectorType> 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);
}; };
struct Node { struct Node {
@ -164,6 +170,7 @@ struct Node {
CoordType initialNormal; CoordType initialNormal;
Vector6d acceleration{0}; Vector6d acceleration{0};
Forces force; Forces force;
Vector6d previousVelocity{0};
Vector6d velocity{0}; Vector6d velocity{0};
double kineticEnergy{0}; double kineticEnergy{0};
Vector6d displacements{0}; Vector6d displacements{0};