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

Average residual forces norm. Alpha angles stored in vectors. Residual forces criterion applicable not only at local kinetic energy maxima

This commit is contained in:
iasonmanolas 2021-07-14 19:17:02 +03:00
parent 5ebf354dcf
commit 8391646d56
4 changed files with 107 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

161
drmsimulationmodel.cpp
View File

@ -222,7 +222,6 @@ void DRMSimulationModel::reset()
mCurrentSimulationStep = 0;
history.clear();
constrainedVertices.clear();
rigidSupports.clear();
pMesh.reset();
plotYValues.clear();
plotHandle.reset();
@ -232,6 +231,8 @@ void DRMSimulationModel::reset()
numOfDampings = 0;
shouldTemporarilyDampForces = false;
externalLoadStep = 1;
isVertexConstrained.clear();
minTotalResidualForcesNorm = std::numeric_limits<double>::max();
}
VectorType DRMSimulationModel::computeDisplacementDifferenceDerivative(
@ -544,13 +545,18 @@ double DRMSimulationModel::computeTheta3(const EdgeType &e, const VertexType &v)
// Use nR as theta3 only for the first star edge
return nR;
}
std::vector<int> incidentElementsIndices(node.incidentElements.size());
for (int iei = 0; iei < incidentElementsIndices.size(); iei++) {
incidentElementsIndices[iei] = pMesh->getIndex(node.incidentElements[iei]);
}
// std::vector<int> incidentElementsIndices(node.incidentElements.size());
// for (int iei = 0; iei < incidentElementsIndices.size(); iei++) {
// incidentElementsIndices[iei] = pMesh->getIndex(node.incidentElements[iei]);
// }
assert(pMesh->getIndex(e) == ei);
assert(node.alphaAngles.contains(ei));
const double alphaAngle = node.alphaAngles.at(elem.ei);
// assert(node.alphaAngles.contains(ei));
const double alphaAngle = std::find_if(node.alphaAngles.begin(),
node.alphaAngles.end(),
[&](const std::pair<EdgeIndex, double> &p) {
return elem.ei == p.first;
})
->second;
const EdgeType &refElem = *node.referenceElement;
const double rotationAngle = nR + alphaAngle;
@ -591,8 +597,7 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e,
{
const Node &node = pMesh->nodes[v];
const VertexIndex &vi = pMesh->nodes[v].vi;
const bool isRigidSupport = rigidSupports.contains(vi);
if (&e == node.referenceElement && !isRigidSupport) {
if (&e == node.referenceElement && !isRigidSupport[vi]) {
if (dui.dofi == DoF::Nr && &dui.v == &v) {
return 1;
} else {
@ -610,7 +615,7 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e,
const VertexPointer &vp_jplus1 = e.cV(1);
double derivativeTheta3_dofi = 0;
if (isRigidSupport) {
if (isRigidSupport[vi]) {
const VectorType &t1Initial = computeT1Vector(pMesh->nodes[vp_j].initialLocation,
pMesh->nodes[vp_jplus1].initialLocation);
VectorType g1 = Cross(t1, t1Initial);
@ -800,9 +805,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
pMesh->EN(), std::vector<std::pair<int, Vector6d>>(4, {-1, Vector6d()}));
// omp_lock_t writelock;
// omp_init_lock(&writelock);
//#ifdef ENABLE_OPENMP
//#pragma omp parallel for schedule(static) num_threads(8)
//#endif
#ifdef ENABLE_OPENMP
#pragma omp parallel for schedule(static) num_threads(5)
#endif
for (int ei = 0; ei < pMesh->EN(); ei++) {
const EdgeType &e = pMesh->edge[ei];
const SimulationMesh::VertexType &ev_j = *e.cV(0);
@ -843,7 +848,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
const size_t vi = edgeNode.vi;
// #pragma omp parallel for schedule(static) num_threads(6)
for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
const bool isDofConstrainedFor_ev = fixedVertices.contains(edgeNode.vi)
const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi]
&& fixedVertices.at(edgeNode.vi).contains(dofi);
if (!isDofConstrainedFor_ev) {
DifferentiateWithRespectTo dui{ev, dofi};
@ -924,7 +929,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
}
if (edgeNode.referenceElement != &e) {
const bool isDofConstrainedFor_refElemOtherVertex
= fixedVertices.contains(refElemOtherVertexNode.vi)
= isVertexConstrained[refElemOtherVertexNode.vi]
&& fixedVertices.at(refElemOtherVertexNode.vi).contains(dofi);
if (!isDofConstrainedFor_refElemOtherVertex) {
DifferentiateWithRespectTo dui{*refElemOtherVertex, dofi};
@ -1012,6 +1017,11 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
}
pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm;
pMesh->totalResidualForcesNorm = totalResidualForcesNorm;
if (mSettings.beVerbose) {
if (minTotalResidualForcesNorm > pMesh->totalResidualForcesNorm) {
minTotalResidualForcesNorm = pMesh->totalResidualForcesNorm;
}
}
pMesh->averageResidualForcesNorm = totalResidualForcesNorm / pMesh->VN();
// pMesh->averageResidualForcesNorm = sumOfResidualForces.norm() / pMesh->VN();
@ -1078,6 +1088,8 @@ void DRMSimulationModel::updateResidualForces()
void DRMSimulationModel::computeRigidSupports()
{
isRigidSupport.clear();
isRigidSupport.resize(pMesh->VN(), false);
for (const VertexType &v : pMesh->vert) {
const VertexIndex vi = pMesh->nodes[v].vi;
const bool isVertexConstrained = constrainedVertices.contains(vi);
@ -1090,7 +1102,7 @@ void DRMSimulationModel::computeRigidSupports()
&& constrainedDoFType.contains(DoF::Ny)
&& constrainedDoFType.contains(DoF::Nr);
if (hasAllDoFTypeConstrained) {
rigidSupports.insert(vi);
isRigidSupport[vi] = true;
}
}
}
@ -1180,8 +1192,9 @@ void DRMSimulationModel::updateNodalMasses()
if (shouldTemporarilyDampForces && mCurrentSimulationStep < untilStep) {
gamma *= 1e6 * (1 - static_cast<double>(mCurrentSimulationStep) / untilStep);
}
if (mCurrentSimulationStep == untilStep && shouldTemporarilyDampForces) {
Dt = mSettings.Dtini * 0.95;
if (mCurrentSimulationStep == static_cast<size_t>(1.2 * untilStep)
&& shouldTemporarilyDampForces) {
Dt = mSettings.Dtini;
}
for (VertexType &v : pMesh->vert) {
const size_t vi = pMesh->getIndex(v);
@ -1328,7 +1341,7 @@ void DRMSimulationModel::updateNodeNr(VertexType &v)
{
const VertexIndex &vi = pMesh->nodes[v].vi;
Node &node = pMesh->nodes[v];
if (!rigidSupports.contains(vi)) {
if (!isRigidSupport[vi]) {
node.nR = node.displacements[5];
} else {
const EdgePointer &refElem = node.referenceElement;
@ -1509,11 +1522,11 @@ void DRMSimulationModel::resetVelocities()
{
for (const VertexType &v : pMesh->vert) {
pMesh->nodes[v].velocity =
// pMesh->nodes[v].acceleration
// * Dt; // NOTE: Do I reset the previous
// velocity?
// reset
// current to 0 or this?
// pMesh->nodes[v].acceleration * Dt
// * 0.5; // NOTE: Do I reset the previous
// // velocity?
// // reset
// // current to 0 or this?
0;
}
updateKineticEnergy();
@ -1616,7 +1629,7 @@ void DRMSimulationModel::updatePositionsOnTheFly(
VectorType newNormal(nx, ny, nz);
v.N() = newNormal;
}
if (!rigidSupports.contains(vi)) {
if (!isRigidSupport[vi]) {
node.nR = node.displacements[5];
} else {
}
@ -1650,15 +1663,18 @@ void DRMSimulationModel::printCurrentState() const
<< std::endl;
std::cout << "Kinetic energy:" << pMesh->currentTotalKineticEnergy << std::endl;
static std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
const auto timePerNodePerIteration = std::chrono::duration_cast<std::chrono::microseconds>(
std::chrono::steady_clock::now() - begin)
.count()
* 1e-6 / (mCurrentSimulationStep * pMesh->VN());
const double timePerNodePerIteration = std::chrono::duration_cast<std::chrono::microseconds>(
std::chrono::steady_clock::now() - begin)
.count()
* 1e-6
/ (static_cast<double>(mCurrentSimulationStep)
* pMesh->VN());
std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" << std::endl;
std::cout << "time(s)/(iterations*node) = " << timePerNodePerIteration << std::endl;
std::cout << "Mov aver deriv norm:" << pMesh->residualForcesMovingAverageDerivativeNorm
<< std::endl;
std::cout << "xi:" << mSettings.xi << std::endl;
std::cout << "Dt:" << Dt << std::endl;
}
void DRMSimulationModel::printDebugInfo() const
@ -1936,7 +1952,7 @@ void DRMSimulationModel::applySolutionGuess(const SimulationResults &solutionGue
node.displacements[5] = nr;
}
// const double nr_asin = q_nr.x()
if (rigidSupports.contains(vi)) {
if (isRigidSupport[vi]) {
const EdgePointer &refElem = node.referenceElement;
const VectorType &refT1 = computeT1Vector(refElem->cP(0), refElem->cP(1));
@ -2046,8 +2062,9 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
}
vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
computeRigidSupports();
isVertexConstrained.resize(pMesh->VN(), false);
for (auto fixedVertex : pJob->constrainedVertices) {
assert(fixedVertex.first < pMesh->VN());
isVertexConstrained[fixedVertex.first] = true;
}
#ifdef POLYSCOPE_DEFINED
@ -2058,7 +2075,8 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
pMesh->getEigenEdges());
polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel,
pMesh->getEigenVertices(),
pMesh->getEigenEdges());
pMesh->getEigenEdges())
->setRadius(0.002);
// registerWorldAxes();
}
#endif
@ -2120,9 +2138,9 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// pJob->save("./PatternOptimizationNonConv");
// Dt = mSettings.Dtini;
}
if (mCurrentSimulationStep == 500 && shouldTemporarilyDampForces) {
Dt = mSettings.Dtini;
}
// if (mCurrentSimulationStep == 500 && shouldTemporarilyDampForces) {
// Dt = mSettings.Dtini;
// }
// while (true) {
updateNormalDerivatives();
updateT1Derivatives();
@ -2208,19 +2226,19 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// = std::chrono::duration_cast<std::chrono::minutes>(t2 - beginTime).count();
// std::cout << "Execution time(min):" << executionTime_min << std::endl;
if (mSettings.useAverage) {
std::cout << "Percentage of target (average):"
std::cout << "Best percentage of target (average):"
<< 100 * mSettings.averageResidualForcesCriterionThreshold
* totalExternalForcesNorm
/ (pMesh->totalResidualForcesNorm / pMesh->VN())
/ (minTotalResidualForcesNorm / pMesh->VN())
<< "%" << std::endl;
}
std::cout << "Percentage of target:"
std::cout << "Best percentage of target:"
<< 100 * mSettings.totalExternalForcesNormPercentageTermination
* totalExternalForcesNorm / pMesh->totalResidualForcesNorm
* totalExternalForcesNorm / minTotalResidualForcesNorm
<< "%" << std::endl;
SimulationResultsReporter::createPlot("Number of Steps",
"Residual Forces mov aver",
history.residualForcesMovingAverage);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver",
// history.residualForcesMovingAverage);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver deriv",
// movingAveragesDerivatives);
@ -2297,8 +2315,18 @@ currentSimulationStep > maxDRMIterations*/
// << std::endl;
// std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm
// << std::endl;
const bool fullfillsResidualForcesNormTerminationCriterion
= !mSettings.useAverage
&& pMesh->totalResidualForcesNorm / totalExternalForcesNorm
< mSettings.totalExternalForcesNormPercentageTermination;
const bool fullfillsAverageResidualForcesNormTerminationCriterion
= mSettings.useAverage
&& (pMesh->totalResidualForcesNorm / pMesh->VN()) / totalExternalForcesNorm
< mSettings.averageResidualForcesCriterionThreshold;
// Residual forces norm convergence
if ((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy
if (((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy
|| fullfillsAverageResidualForcesNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion)
// && mCurrentSimulationStep > movingAverageSampleSize
&& (pJob->nodalForcedDisplacements.empty()
|| mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps))
@ -2314,13 +2342,6 @@ mesh->currentTotalPotentialEnergykN*/
&& pMesh->currentTotalTranslationalKineticEnergy
< mSettings.totalTranslationalKineticEnergyThreshold
&& mCurrentSimulationStep > 20 && numOfDampings > 0;
const bool fullfillsResidualForcesNormTerminationCriterion
= pMesh->totalResidualForcesNorm / totalExternalForcesNorm
< mSettings.totalExternalForcesNormPercentageTermination;
const bool fullfillsAverageResidualForcesNormTerminationCriterion
= mSettings.useAverage
&& (pMesh->totalResidualForcesNorm / pMesh->VN()) / totalExternalForcesNorm
< mSettings.averageResidualForcesCriterionThreshold;
const bool fullfillsMovingAverageNormTerminationCriterion
= pMesh->residualForcesMovingAverage
< mSettings.residualForcesMovingAverageNormThreshold;
@ -2341,7 +2362,7 @@ mesh->currentTotalPotentialEnergykN*/
std::cout << "Converged using residual forces norm threshold criterion"
<< std::endl;
} else if (fullfillsKineticEnergyTerminationCriterion) {
std::cout << "Warning: The kinetic energy of the system was "
std::cout << "The kinetic energy of the system was "
" used as a convergence criterion"
<< std::endl;
} else if (fullfillsMovingAverageNormTerminationCriterion) {
@ -2377,25 +2398,25 @@ mesh->currentTotalPotentialEnergykN*/
// }
}
// 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 * 0.5 * 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;
@ -2408,11 +2429,11 @@ mesh->currentTotalPotentialEnergykN*/
// + xi_init - triggerPercentage * xi_min)
// / (1 - triggerPercentage);
// }
resetVelocities();
Dt *= mSettings.xi;
// if (mSettings.isDebugMode) {
// std::cout << Dt << std::endl;
// }
resetVelocities();
++numOfDampings;
if (mSettings.shouldCreatePlots) {
history.redMarks.push_back(mCurrentSimulationStep);

11
drmsimulationmodel.hpp
View File

@ -45,7 +45,7 @@ public:
double totalResidualForcesNormThreshold{1e-3};
double totalExternalForcesNormPercentageTermination{1e-3};
bool useAverage{false};
double averageResidualForcesCriterionThreshold{1e-3};
double averageResidualForcesCriterionThreshold{1e-5};
Settings() {}
};
@ -61,6 +61,9 @@ private:
std::vector<double> plotYValues;
size_t numOfDampings{0};
int externalLoadStep{1};
std::vector<bool> isVertexConstrained;
std::vector<bool> isRigidSupport;
double minTotalResidualForcesNorm{std::numeric_limits<double>::max()};
const std::string meshPolyscopeLabel{"Simulation mesh"};
std::unique_ptr<SimulationMesh> pMesh;
@ -69,7 +72,6 @@ private:
// Eigen::Tensor<double, 4> theta3Derivatives;
// std::unordered_map<MyKeyType, double, key_hash> theta3Derivatives;
bool shouldApplyInitialDistortion = false;
std::unordered_set<VertexIndex> rigidSupports;
void reset();
void updateNodalInternalForces(
@ -156,10 +158,9 @@ private:
void updateElementalFrames();
VectorType computeDerivativeOfR(const EdgeType &e,
const DifferentiateWithRespectTo &dui) const;
VectorType computeDerivativeOfR(const EdgeType &e, const DifferentiateWithRespectTo &dui) const;
bool isRigidSupport(const VertexType &v) const;
// bool isRigidSupport(const VertexType &v) const;
static double computeDerivativeOfNorm(const VectorType &x,
const VectorType &derivativeOfX);

8
simulationmesh.cpp
View File

@ -140,6 +140,7 @@ void SimulationMesh::initializeNodes() {
const EdgeType &referenceElement = *node.referenceElement;
const VectorType t01 = computeT1Vector(referenceElement.cP(0), referenceElement.cP(1));
const VectorType f01 = (t01 - (v.cN() * (t01.dot(v.cN())))).Normalize();
node.alphaAngles.reserve(incidentElements.size());
for (const VCGEdgeMesh::EdgePointer &ep : nodes[v].incidentElements) {
assert(referenceElement.cV(0) == ep->cV(0) || referenceElement.cV(0) == ep->cV(1)
@ -148,7 +149,7 @@ void SimulationMesh::initializeNodes() {
const VectorType f1 = t1 - (v.cN() * (t1.dot(v.cN()))).Normalize();
const EdgeIndex ei = getIndex(ep);
const double alphaAngle = computeAngle(f01, f1, v.cN());
node.alphaAngles[ei] = alphaAngle;
node.alphaAngles.emplace_back(std::make_pair(ei, alphaAngle));
}
}
}
@ -182,7 +183,8 @@ void SimulationMesh::reset() {
const EdgeType &referenceElement = *getReferenceElement(v);
const VectorType t01 = computeT1Vector(referenceElement.cP(0), referenceElement.cP(1));
const VectorType f01 = (t01 - (v.cN() * (t01.dot(v.cN())))).Normalize();
node.alphaAngles.clear();
node.alphaAngles.reserve(node.incidentElements.size());
for (const VCGEdgeMesh::EdgePointer &ep : nodes[v].incidentElements) {
assert(referenceElement.cV(0) == ep->cV(0) || referenceElement.cV(0) == ep->cV(1)
|| referenceElement.cV(1) == ep->cV(0) || referenceElement.cV(1) == ep->cV(1));
@ -190,7 +192,7 @@ void SimulationMesh::reset() {
const VectorType f1 = t1 - (v.cN() * (t1.dot(v.cN()))).Normalize();
const EdgeIndex ei = getIndex(ep);
const double alphaAngle = computeAngle(f01, f1, v.cN());
node.alphaAngles[ei] = alphaAngle;
node.alphaAngles.emplace_back(std::make_pair(ei, alphaAngle));
}
}
}

9
simulationmesh.hpp
View File

@ -157,10 +157,11 @@ struct Node {
double kineticEnergy{0};
Vector6d displacements{0};
double nR{0};
std::unordered_map<EdgeIndex, double>
alphaAngles; // contains the initial angles between the first star element
// incident to this node and the other elements of the star
// has size equal to the valence of the vertex
// std::unordered_map<EdgeIndex, double>
// alphaAngles; // contains the initial angles between the first star element
// // incident to this node and the other elements of the star
// // has size equal to the valence of the vertex
std::vector<std::pair<EdgeIndex, double>> alphaAngles;
std::vector<VCGEdgeMesh::EdgePointer> incidentElements;
std::vector<VectorType> derivativeOfNormal;