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

Merge branch 'drmGuess'

This commit is contained in:
iasonmanolas 2021-06-24 10:03:23 +03:00
parent e5ba83c0c8 3e3d32399b
commit 9b6aba7bac
15 changed files with 829 additions and 412 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

730
drmsimulationmodel.cpp
View File

@ -45,10 +45,10 @@ void DRMSimulationModel::runUnitTests()
settings.Dtini = 0.1; settings.Dtini = 0.1;
settings.xi = 0.9969; settings.xi = 0.9969;
settings.maxDRMIterations = 0; settings.maxDRMIterations = 0;
formFinder.totalResidualForcesNormThreshold = 1; formFinder.mSettings.totalResidualForcesNormThreshold = 1;
settings.shouldDraw = false; settings.shouldDraw = false;
settings.beVerbose = false; settings.beVerbose = 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);
simpleBeam_simulationResults.save(); simpleBeam_simulationResults.save();
@ -227,11 +227,10 @@ void DRMSimulationModel::reset()
plotYValues.clear(); plotYValues.clear();
plotHandle.reset(); plotHandle.reset();
checkedForMaximumMoment = false; checkedForMaximumMoment = false;
mSettings.shouldUseTranslationalKineticEnergyThreshold = false;
externalMomentsNorm = 0; externalMomentsNorm = 0;
mSettings.drawingStep = 1;
Dt = mSettings.Dtini; Dt = mSettings.Dtini;
numOfDampings = 0; numOfDampings = 0;
shouldTemporarilyDampForces = false;
} }
VectorType DRMSimulationModel::computeDisplacementDifferenceDerivative( VectorType DRMSimulationModel::computeDisplacementDifferenceDerivative(
@ -539,6 +538,7 @@ double DRMSimulationModel::computeTheta3(const EdgeType &e, const VertexType &v)
// Should be refactored in the future // Should be refactored in the future
double theta3; double theta3;
const bool shouldBreak = mCurrentSimulationStep == 12970;
if (&e == node.referenceElement) { if (&e == node.referenceElement) {
// Use nR as theta3 only for the first star edge // Use nR as theta3 only for the first star edge
return nR; return nR;
@ -794,7 +794,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
// omp_lock_t writelock; // omp_lock_t writelock;
// omp_init_lock(&writelock); // omp_init_lock(&writelock);
#ifdef ENABLE_OPENMP #ifdef ENABLE_OPENMP
#pragma omp parallel for schedule(static) num_threads(5) #pragma omp parallel for schedule(static) num_threads(8)
#endif #endif
for (int ei = 0; ei < pMesh->EN(); ei++) { for (int ei = 0; ei < pMesh->EN(); ei++) {
const EdgeType &e = pMesh->edge[ei]; const EdgeType &e = pMesh->edge[ei];
@ -808,6 +808,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
const double theta1_jplus1 = ef.t1.dot(t3CrossN_jplus1); const double theta1_jplus1 = ef.t1.dot(t3CrossN_jplus1);
const double theta2_j = ef.t2.dot(t3CrossN_j); const double theta2_j = ef.t2.dot(t3CrossN_j);
const double theta2_jplus1 = ef.t2.dot(t3CrossN_jplus1); const double theta2_jplus1 = ef.t2.dot(t3CrossN_jplus1);
const bool shouldBreak = mCurrentSimulationStep == 12970;
const double theta3_j = computeTheta3(e, ev_j); const double theta3_j = computeTheta3(e, ev_j);
const double theta3_jplus1 = computeTheta3(e, ev_jplus1); const double theta3_jplus1 = computeTheta3(e, ev_jplus1);
// element.rotationalDisplacements_j.theta1 = theta1_j; // element.rotationalDisplacements_j.theta1 = theta1_j;
@ -837,7 +838,6 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
const bool isDofConstrainedFor_ev = fixedVertices.contains(edgeNode.vi) const bool isDofConstrainedFor_ev = fixedVertices.contains(edgeNode.vi)
&& fixedVertices.at(edgeNode.vi).contains(dofi); && fixedVertices.at(edgeNode.vi).contains(dofi);
const bool shouldBreak = edgeNode.vi == 0 && dofi == 5;
if (!isDofConstrainedFor_ev) { if (!isDofConstrainedFor_ev) {
DifferentiateWithRespectTo dui{ev, dofi}; DifferentiateWithRespectTo dui{ev, dofi};
// Axial force computation // Axial force computation
@ -976,8 +976,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
for (size_t vi = 0; vi < pMesh->VN(); vi++) { for (size_t vi = 0; vi < pMesh->VN(); vi++) {
Node::Forces &force = pMesh->nodes[vi].force; Node::Forces &force = pMesh->nodes[vi].force;
const Vector6d &nodeResidualForce = force.residual; const Vector6d &nodeResidualForce = force.residual;
// sumOfResidualForces = sumOfResidualForces + nodeResidualForce; sumOfResidualForces = sumOfResidualForces + nodeResidualForce;
const double residualForceNorm = nodeResidualForce.norm(); // const double residualForceNorm = nodeResidualForce.norm();
const double residualForceNorm = nodeResidualForce.getTranslation().norm();
const bool shouldTrigger = std::isnan(residualForceNorm); const bool shouldTrigger = std::isnan(residualForceNorm);
totalResidualForcesNorm += residualForceNorm; totalResidualForcesNorm += residualForceNorm;
} }
@ -1145,7 +1146,11 @@ DRMSimulationModel::DRMSimulationModel() {}
void DRMSimulationModel::updateNodalMasses() void DRMSimulationModel::updateNodalMasses()
{ {
const double gamma = 0.8; double gamma = mSettings.gamma;
const size_t untilStep = 500;
if (shouldTemporarilyDampForces && mCurrentSimulationStep < untilStep) {
gamma *= 1e6 * (1 - static_cast<double>(mCurrentSimulationStep) / untilStep);
}
for (VertexType &v : pMesh->vert) { for (VertexType &v : pMesh->vert) {
const size_t vi = pMesh->getIndex(v); const size_t vi = pMesh->getIndex(v);
double translationalSumSk = 0; double translationalSumSk = 0;
@ -1208,6 +1213,14 @@ void DRMSimulationModel::updateNodalAccelerations()
node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalI2; node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalI2;
} }
} }
// if (vi == 10) {
// std::cout << "Acceleration:" << node.acceleration[0] << " " << node.acceleration[1]
// << " " << node.acceleration[2] << std::endl;
// }
// if (shouldTemporarilyDampForces && mCurrentSimulationStep < 700) {
// node.acceleration = node.acceleration * 1e-2;
// }
} }
} }
@ -1217,25 +1230,35 @@ void DRMSimulationModel::updateNodalVelocities()
const VertexIndex vi = pMesh->getIndex(v); const VertexIndex vi = pMesh->getIndex(v);
Node &node = pMesh->nodes[v]; Node &node = pMesh->nodes[v];
node.velocity = node.velocity + node.acceleration * Dt; node.velocity = node.velocity + node.acceleration * Dt;
// if (vi == 10) {
// std::cout << "Velocity:" << node.velocity[0] << " " << node.velocity[1] << " "
// << node.velocity[2] << std::endl;
// }
} }
updateKineticEnergy(); updateKineticEnergy();
} }
void DRMSimulationModel::updateNodalDisplacements() void DRMSimulationModel::updateNodalDisplacements()
{ {
const bool shouldCapDisplacements = mSettings.displacementCap.has_value();
for (VertexType &v : pMesh->vert) { for (VertexType &v : pMesh->vert) {
Node &node = pMesh->nodes[v]; Node &node = pMesh->nodes[v];
node.displacements = node.displacements + node.velocity * Dt; Vector6d stepDisplacement = node.velocity * Dt;
// if (mSettings.beVerbose) { if (shouldCapDisplacements
// std::cout << "Node " << node.vi << ":" << endl; && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) {
// std::cout << node.displacements[0] << " " << node.displacements[0] stepDisplacement = stepDisplacement
// << " " * (*mSettings.displacementCap
// << node.displacements[1] << " " << node.displacements[2] / stepDisplacement.getTranslation().norm());
// << " " std::cout << "Displacement capped" << std::endl;
// << node.displacements[3] << " " << node.displacements[4] }
// << " " node.displacements = node.displacements + stepDisplacement;
// << node.displacements[5] << std::endl; // if (mSettings.isDebugMode && mSettings.beVerbose && pMesh->getIndex(v) == 40) {
// } // std::cout << "Node " << node.vi << ":" << endl;
// std::cout << node.displacements[0] << " " << node.displacements[0] << " "
// << node.displacements[1] << " " << node.displacements[2] << " "
// << node.displacements[3] << " " << node.displacements[4] << " "
// << node.displacements[5] << std::endl;
// }
} }
} }
@ -1263,6 +1286,23 @@ void DRMSimulationModel::updateNodePosition(
} }
} }
void DRMSimulationModel::updateNodeNr(VertexType &v)
{
const VertexIndex &vi = pMesh->nodes[v].vi;
Node &node = pMesh->nodes[v];
if (!rigidSupports.contains(vi)) {
node.nR = node.displacements[5];
} else {
const EdgePointer &refElem = node.referenceElement;
const VectorType &refT1 = pMesh->elements[refElem].frame.t1;
const VectorType &t1Initial = computeT1Vector(pMesh->nodes[refElem->cV(0)].initialLocation,
pMesh->nodes[refElem->cV(1)].initialLocation);
VectorType g1 = Cross(refT1, t1Initial);
node.nR = g1.dot(v.cN());
}
}
void DRMSimulationModel::updateNodeNormal( void DRMSimulationModel::updateNodeNormal(
VertexType &v, const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices) VertexType &v, const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices)
{ {
@ -1319,19 +1359,6 @@ void DRMSimulationModel::updateNodeNormal(
} }
// if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(DoF::Nr)) { // if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(DoF::Nr)) {
if (!rigidSupports.contains(vi)) {
node.nR = node.displacements[5];
} else {
const EdgePointer &refElem = node.referenceElement;
const VectorType &refT1 = pMesh->elements[refElem].frame.t1;
const VectorType &t1Initial = computeT1Vector(pMesh->nodes[refElem->cV(0)].initialLocation,
pMesh->nodes[refElem->cV(1)].initialLocation);
VectorType g1 = Cross(refT1, t1Initial);
node.nR = g1.dot(v.cN());
int i = 0;
i++;
}
// if (vi == 1) { // if (vi == 1) {
// std::cout << "AFTER:" << mesh->vert[1].N()[0] << " " << // std::cout << "AFTER:" << mesh->vert[1].N()[0] << " " <<
// mesh->vert[1].N()[1] // mesh->vert[1].N()[1]
@ -1345,6 +1372,7 @@ void DRMSimulationModel::applyDisplacements(
for (VertexType &v : pMesh->vert) { for (VertexType &v : pMesh->vert) {
updateNodePosition(v, fixedVertices); updateNodePosition(v, fixedVertices);
updateNodeNormal(v, fixedVertices); updateNodeNormal(v, fixedVertices);
updateNodeNr(v);
} }
updateElementalFrames(); updateElementalFrames();
if (mSettings.shouldDraw) { if (mSettings.shouldDraw) {
@ -1427,8 +1455,8 @@ void DRMSimulationModel::updateKineticEnergy()
const double nodeRotationalKineticEnergy const double nodeRotationalKineticEnergy
= 0.5 = 0.5
* (node.mass.rotationalJ * pow(node.velocity[3], 2) * (node.mass.rotationalJ * pow(node.velocity[3], 2)
+ +node.mass.rotationalI3 * pow(node.velocity[4], 2) + node.mass.rotationalI3 * pow(node.velocity[4], 2)
+ +node.mass.rotationalI2 * pow(node.velocity[5], 2)); + node.mass.rotationalI2 * pow(node.velocity[5], 2));
node.kineticEnergy += nodeTranslationalKineticEnergy /*+ nodeRotationalKineticEnergy*/; node.kineticEnergy += nodeTranslationalKineticEnergy /*+ nodeRotationalKineticEnergy*/;
assert(node.kineticEnergy < 1e15); assert(node.kineticEnergy < 1e15);
@ -1442,12 +1470,13 @@ void DRMSimulationModel::updateKineticEnergy()
void DRMSimulationModel::resetVelocities() void DRMSimulationModel::resetVelocities()
{ {
for (const VertexType &v : pMesh->vert) { for (const VertexType &v : pMesh->vert) {
pMesh->nodes[v].velocity = 0; pMesh->nodes[v].velocity =
// mesh->nodes[v].force.residual * 0.5 * Dt / // pMesh->nodes[v].acceleration
// mesh->nodes[v].mass; // NOTE: Do I reset the previous // * Dt; // NOTE: Do I reset the previous
// velocity? // velocity?
// reset // reset
// current to 0 or this? // current to 0 or this?
0;
} }
updateKineticEnergy(); updateKineticEnergy();
} }
@ -1584,6 +1613,9 @@ void DRMSimulationModel::printCurrentState() const
* 1e-6 / (mCurrentSimulationStep * pMesh->VN()); * 1e-6 / (mCurrentSimulationStep * pMesh->VN());
std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" << std::endl; std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" << std::endl;
std::cout << "time(s)/(iterations*node) = " << timePerNodePerIteration << 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;
} }
void DRMSimulationModel::printDebugInfo() const void DRMSimulationModel::printDebugInfo() const
@ -1660,12 +1692,13 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
accelerationX[pMesh->getIndex(v)] = pMesh->nodes[v].acceleration[0]; accelerationX[pMesh->getIndex(v)] = pMesh->nodes[v].acceleration[0];
} }
meshPolyscopeHandle->addNodeScalarQuantity("Kinetic Energy", kineticEnergies)->setEnabled(false); meshPolyscopeHandle->addNodeScalarQuantity("Kinetic Energy", kineticEnergies)->setEnabled(false);
meshPolyscopeHandle->addNodeVectorQuantity("Node normals", nodeNormals)->setEnabled(true); meshPolyscopeHandle->addNodeVectorQuantity("Node normals", nodeNormals)->setEnabled(false);
meshPolyscopeHandle->addNodeVectorQuantity("Internal force", internalForces)->setEnabled(false); meshPolyscopeHandle->addNodeVectorQuantity("Internal force", internalForces)->setEnabled(false);
meshPolyscopeHandle->addNodeVectorQuantity("External force", externalForces)->setEnabled(false); meshPolyscopeHandle->addNodeVectorQuantity("External force", externalForces)->setEnabled(false);
meshPolyscopeHandle->addNodeVectorQuantity("External Moments", externalMoments)->setEnabled(true); meshPolyscopeHandle->addNodeVectorQuantity("External Moments", externalMoments)->setEnabled(true);
meshPolyscopeHandle->addNodeScalarQuantity("Internal force norm", internalForcesNorm) meshPolyscopeHandle->addNodeScalarQuantity("Internal force norm", internalForcesNorm)
->setEnabled(true); ->setEnabled(true);
meshPolyscopeHandle->setRadius(0.002);
// polyscope::getCurveNetwork(meshPolyscopeLabel) // polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeVectorQuantity("Internal Axial force", internalAxialForces) // ->addNodeVectorQuantity("Internal Axial force", internalAxialForces)
// ->setEnabled(false); // ->setEnabled(false);
@ -1735,8 +1768,9 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
// instead of full width. Must have // instead of full width. Must have
// matching PopItemWidth() below. // matching PopItemWidth() below.
ImGui::InputInt("Simulation drawing step", ImGui::InputInt("Simulation debug step",
&mSettings.drawingStep); // set a int variable &mSettings.debugModeStep); // set a int variable
mSettings.drawingStep = mSettings.debugModeStep;
ImGui::Checkbox("Enable drawing", ImGui::Checkbox("Enable drawing",
&mSettings.shouldDraw); // set a int variable &mSettings.shouldDraw); // set a int variable
ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep); ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep);
@ -1766,30 +1800,183 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
} }
#endif #endif
void DRMSimulationModel::applySolutionGuess(const SimulationResults &solutionGuess,
const std::shared_ptr<SimulationJob> &pJob)
{
assert(solutionGuess.displacements.size() == pMesh->VN()
&& solutionGuess.rotationalDisplacementQuaternion.size() == pMesh->VN());
for (size_t vi = 0; vi < pMesh->VN(); vi++) {
Node &node = pMesh->nodes[vi];
Eigen::Vector3d translationalDisplacements(solutionGuess.displacements[vi].getTranslation());
if (!pJob->constrainedVertices.contains(vi)
|| !pJob->constrainedVertices.at(vi).contains(0)) {
node.displacements[0] = translationalDisplacements[0];
}
if (!pJob->constrainedVertices.contains(vi)
|| !pJob->constrainedVertices.at(vi).contains(1)) {
node.displacements[1] = translationalDisplacements[1];
}
if (!pJob->constrainedVertices.contains(vi)
|| !pJob->constrainedVertices.at(vi).contains(2)) {
node.displacements[2] = translationalDisplacements[2];
}
updateNodePosition(pMesh->vert[vi], pJob->constrainedVertices);
}
for (size_t vi = 0; vi < pMesh->VN(); vi++) {
Node &node = pMesh->nodes[vi];
const Eigen::Vector3d nInitial_eigen = node.initialNormal.ToEigenVector<Eigen::Vector3d>();
Eigen::Quaternion<double> q;
Eigen::Vector3d eulerAngles = solutionGuess.displacements[vi].getRotation();
q = Eigen::AngleAxisd(eulerAngles[0], Eigen::Vector3d::UnitX())
* Eigen::AngleAxisd(eulerAngles[1], Eigen::Vector3d::UnitY())
* Eigen::AngleAxisd(eulerAngles[2], Eigen::Vector3d::UnitZ());
Eigen::Vector3d nDeformed_eigen = (q * nInitial_eigen) /*.normalized()*/;
nDeformed_eigen.normalized();
// Eigen::Vector3d n_groundOfTruth = solutionGuess.debug_q_normal[vi] * nInitial_eigen;
// n_groundOfTruth.normalized();
if (!pJob->constrainedVertices.contains(vi)
|| !pJob->constrainedVertices.at(vi).contains(3)) {
node.displacements[3] = (nDeformed_eigen - nInitial_eigen)[0];
}
if (!pJob->constrainedVertices.contains(vi)
|| !pJob->constrainedVertices.at(vi).contains(4)) {
node.displacements[4] = (nDeformed_eigen - nInitial_eigen)[1];
}
updateNodeNormal(pMesh->vert[vi], pJob->constrainedVertices);
// const auto debug_rightNy = solutionGuess.debug_drmDisplacements[vi][4];
Eigen::Vector3d referenceT1_deformed = computeT1Vector(node.referenceElement->cP(0),
node.referenceElement->cP(1))
.ToEigenVector<Eigen::Vector3d>();
const Eigen::Vector3d referenceF1_deformed
= (referenceT1_deformed - (nInitial_eigen * (referenceT1_deformed.dot(nInitial_eigen))))
.normalized();
const Eigen::Vector3d referenceT1_initial
= computeT1Vector(pMesh->nodes[node.referenceElement->cV(0)].initialLocation,
pMesh->nodes[node.referenceElement->cV(1)].initialLocation)
.ToEigenVector<Eigen::Vector3d>();
// const VectorType &n_initial = node.initialNormal;
const Eigen::Vector3d referenceF1_initial
= (referenceT1_initial - (nInitial_eigen * (referenceT1_initial.dot(nInitial_eigen))))
.normalized();
Eigen::Quaternion<double> q_f1; //nr is with respect to f1
q_f1.setFromTwoVectors(referenceF1_initial, referenceF1_deformed);
Eigen::Quaternion<double> q_normal;
q_normal.setFromTwoVectors(nInitial_eigen, nDeformed_eigen);
Eigen::Quaternion<double> q_nr = q_f1.inverse() * q_normal.inverse() * q;
const double nr_0To2pi_pos = 2 * std::acos(q_nr.w());
// const double nr_0To2pi_groundOfTruth = 2 * std::acos(solutionGuess.debug_q_nr[vi].w());
const double nr_0To2pi = nr_0To2pi_pos <= M_PI ? nr_0To2pi_pos : (nr_0To2pi_pos - 2 * M_PI);
Eigen::Vector3d deformedNormal_debug(q_nr.x() * sin(nr_0To2pi_pos / 2),
q_nr.y() * sin(nr_0To2pi_pos / 2),
q_nr.z() * sin(nr_0To2pi_pos / 2));
/*deformedNormal_debug =*/deformedNormal_debug.normalize();
const double nr = deformedNormal_debug.dot(nDeformed_eigen) > 0 ? nr_0To2pi : -nr_0To2pi;
if (!pJob->constrainedVertices.contains(vi)
|| !pJob->constrainedVertices.at(vi).contains(5)) {
node.displacements[5] = nr;
}
// const double nr_asin = q_nr.x()
if (rigidSupports.contains(vi)) {
const EdgePointer &refElem = node.referenceElement;
const VectorType &refT1 = computeT1Vector(refElem->cP(0), refElem->cP(1));
const VectorType &t1Initial
= computeT1Vector(pMesh->nodes[refElem->cV(0)].initialLocation,
pMesh->nodes[refElem->cV(1)].initialLocation);
VectorType g1 = Cross(refT1, t1Initial);
node.nR = g1.dot(pMesh->vert[vi].cN());
} else {
node.displacements[5] = nr;
}
}
updateElementalFrames();
applyDisplacements(constrainedVertices);
if (!pJob->nodalForcedDisplacements.empty()) {
applyForcedDisplacements(
pJob->nodalForcedDisplacements);
}
updateElementalLengths();
// // registerWorldAxes();
// Eigen::MatrixX3d framesX(pMesh->VN(), 3);
// Eigen::MatrixX3d framesY(pMesh->VN(), 3);
// Eigen::MatrixX3d framesZ(pMesh->VN(), 3);
// for (VertexIndex vi = 0; vi < pMesh->VN(); vi++) {
// Node &node = pMesh->nodes[vi];
// Eigen::Vector3d translationalDisplacements(solutionGuess.displacements[vi].getTranslation());
// node.displacements[0] = translationalDisplacements[0];
// node.displacements[1] = translationalDisplacements[1];
// node.displacements[2] = translationalDisplacements[2];
// Eigen::Quaternion<double> q;
// Eigen::Vector3d eulerAngles = solutionGuess.displacements[vi].getRotation();
// q = Eigen::AngleAxisd(eulerAngles[0], Eigen::Vector3d::UnitX())
// * Eigen::AngleAxisd(eulerAngles[1], Eigen::Vector3d::UnitY())
// * Eigen::AngleAxisd(eulerAngles[2], Eigen::Vector3d::UnitZ());
// auto deformedNormal = q * Eigen::Vector3d(0, 0, 1);
// auto deformedFrameY = q * Eigen::Vector3d(0, 1, 0);
// auto deformedFrameX = q * Eigen::Vector3d(1, 0, 0);
// framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0], deformedFrameX[1], deformedFrameX[2]);
// framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0], deformedFrameY[1], deformedFrameY[2]);
// framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0], deformedNormal[1], deformedNormal[2]);
// }
// polyscope::CurveNetwork *meshPolyscopeHandle = polyscope::getCurveNetwork(meshPolyscopeLabel);
// meshPolyscopeHandle->updateNodePositions(pMesh->getEigenVertices());
// //if(showFramesOn.contains(vi)){
// auto polyscopeHandle_frameX = meshPolyscopeHandle->addNodeVectorQuantity("FrameX", framesX);
// polyscopeHandle_frameX->setVectorLengthScale(0.01);
// polyscopeHandle_frameX->setVectorRadius(0.01);
// polyscopeHandle_frameX->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
// auto polyscopeHandle_frameY = meshPolyscopeHandle->addNodeVectorQuantity("FrameY", framesY);
// polyscopeHandle_frameY->setVectorLengthScale(0.01);
// polyscopeHandle_frameY->setVectorRadius(0.01);
// polyscopeHandle_frameY->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
// auto polyscopeHandle_frameZ = meshPolyscopeHandle->addNodeVectorQuantity("FrameZ", framesZ);
// polyscopeHandle_frameZ->setVectorLengthScale(0.01);
// polyscopeHandle_frameZ->setVectorRadius(0.01);
// polyscopeHandle_frameZ->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
// //}
// polyscope::show();
}
SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob, SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings) const Settings &settings,
const SimulationResults &solutionGuess)
{ {
assert(pJob->pMesh.operator bool()); assert(pJob->pMesh.operator bool());
auto t1 = std::chrono::high_resolution_clock::now(); auto t1 = std::chrono::high_resolution_clock::now();
reset();
mSettings = settings; mSettings = settings;
double totalExternalForcesNorm = 0; reset();
for (const std::pair<VertexIndex, Vector6d> &externalForce : pJob->nodalExternalForces) { // double totalExternalForcesNorm = 0;
totalExternalForcesNorm += externalForce.second.norm(); // for (const std::pair<VertexIndex, Vector6d> &externalForce : pJob->nodalExternalForces) {
} // totalExternalForcesNorm += externalForce.second.norm();
// }
if (!pJob->nodalExternalForces.empty()) { // if (!pJob->nodalExternalForces.empty()) {
totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination // // totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination
* totalExternalForcesNorm; // // * totalExternalForcesNorm;
} else { // } else {
totalResidualForcesNormThreshold = 1; // totalResidualForcesNormThreshold = 1;
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:" << totalResidualForcesNormThreshold std::cout << "totalResidualForcesNormThreshold:"
<< std::endl; << mSettings.totalResidualForcesNormThreshold << std::endl;
} }
history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
// if (!pJob->nodalExternalForces.empty()) { // if (!pJob->nodalExternalForces.empty()) {
// double externalForcesNorm = 0; // double externalForcesNorm = 0;
@ -1835,141 +2022,170 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
} }
#endif #endif
updateElementalFrames(); if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) {
updateNodalMasses(); shouldApplyInitialDistortion = true;
if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) { }
shouldApplyInitialDistortion = true; updateElementalFrames();
} if (!solutionGuess.displacements.empty()) {
updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices); //#ifdef POLYSCOPE_DEFINED
const size_t movingAverageSampleSize = 50; // if (mSettings.shouldDraw || mSettings.isDebugMode) {
std::queue<double> residualForcesMovingAverageHistorySample; // solutionGuess.registerForDrawing();
// polyscope::show();
// solutionGuess.unregister();
// }
//#endif
double residualForcesMovingAverageDerivativeMax = 0; applySolutionGuess(solutionGuess, pJob);
while (settings.maxDRMIterations == 0 || mCurrentSimulationStep < settings.maxDRMIterations) { shouldTemporarilyDampForces = true;
// while (true) { }
updateNormalDerivatives();
updateT1Derivatives();
updateRDerivatives();
updateT2Derivatives();
updateT3Derivatives();
updateResidualForcesOnTheFly(constrainedVertices);
// TODO: write parallel function for updating positions updateNodalMasses();
// TODO: Make a single function out of updateResidualForcesOnTheFly updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
// updatePositionsOnTheFly const size_t movingAverageSampleSize = 200;
// updatePositionsOnTheFly(constrainedVertices); std::queue<double> residualForcesMovingAverageHistorySample;
updateNodalMasses(); std::vector<double> percentageOfConvergence;
updateNodalAccelerations();
updateNodalVelocities();
updateNodalDisplacements();
applyDisplacements(constrainedVertices);
if (!pJob->nodalForcedDisplacements.empty()) {
applyForcedDisplacements(
pJob->nodalForcedDisplacements); double residualForcesMovingAverageDerivativeMax = 0;
} while (settings.maxDRMIterations == 0 || mCurrentSimulationStep < settings.maxDRMIterations) {
// if (!pJob->nodalForcedNormals.empty()) { // while (true) {
// applyForcedNormals(pJob->nodalForcedNormals); updateNormalDerivatives();
// } updateT1Derivatives();
updateElementalLengths(); updateRDerivatives();
mCurrentSimulationStep++; updateT2Derivatives();
updateT3Derivatives();
const bool shouldBreak = mCurrentSimulationStep == 12970;
updateResidualForcesOnTheFly(constrainedVertices);
double sumOfDisplacementsNorm = 0; // TODO: write parallel function for updating positions
for (size_t vi = 0; vi < pMesh->VN(); vi++) { // TODO: Make a single function out of updateResidualForcesOnTheFly
sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm(); // updatePositionsOnTheFly
} // updatePositionsOnTheFly(constrainedVertices);
sumOfDisplacementsNorm /= pMesh->bbox.Diag(); updateNodalMasses();
pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm; updateNodalAccelerations();
// pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm updateNodalVelocities();
// + mCurrentSimulationStep updateNodalDisplacements();
// * pMesh->residualForcesMovingAverage) applyDisplacements(constrainedVertices);
// / (1 + mCurrentSimulationStep); if (!pJob->nodalForcedDisplacements.empty()) {
pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm applyForcedDisplacements(
/ movingAverageSampleSize;
residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage);
if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) {
const double firstElementValue = residualForcesMovingAverageHistorySample.front();
pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize;
residualForcesMovingAverageHistorySample.pop();
pMesh->residualForcesMovingAverageDerivativeNorm pJob->nodalForcedDisplacements);
= std::abs((residualForcesMovingAverageHistorySample.back() }
- residualForcesMovingAverageHistorySample.front())) // if (!pJob->nodalForcedNormals.empty()) {
/ (movingAverageSampleSize); // applyForcedNormals(pJob->nodalForcedNormals);
residualForcesMovingAverageDerivativeMax // }
= std::max(pMesh->residualForcesMovingAverageDerivativeNorm, updateElementalLengths();
residualForcesMovingAverageDerivativeMax); mCurrentSimulationStep++;
pMesh->residualForcesMovingAverageDerivativeNorm
/= residualForcesMovingAverageDerivativeMax;
// std::cout << "Normalized derivative:"
// << residualForcesMovingAverageDerivativeNorm
// << std::endl;
}
// pMesh->previousTotalPotentialEnergykN = double sumOfDisplacementsNorm = 0;
// pMesh->currentTotalPotentialEnergykN; for (size_t vi = 0; vi < pMesh->VN(); vi++) {
// pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm();
// timePerNodePerIterationHistor.push_back(timePerNodePerIteration); }
if (mSettings.beVerbose && mSettings.isDebugMode sumOfDisplacementsNorm /= pMesh->bbox.Diag();
&& mCurrentSimulationStep % mSettings.debugModeStep == 0) { pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm;
printCurrentState(); // pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm
// SimulationResultsReporter::createPlot("Number of Steps", // + mCurrentSimulationStep
// "Residual Forces mov aver", // * pMesh->residualForcesMovingAverage)
// movingAverages); // / (1 + mCurrentSimulationStep);
// SimulationResultsReporter::createPlot("Number of Steps", pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm
// "Residual Forces mov aver deriv", / movingAverageSampleSize;
// movingAveragesDerivatives); residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage);
// draw(); if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) {
std::cout << "Mov aver deriv:" << pMesh->residualForcesMovingAverageDerivativeNorm const double firstElementValue = residualForcesMovingAverageHistorySample.front();
<< std::endl; pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize;
// SimulationResulnodalForcedDisplacementstsReporter::createPlot("Number of Steps", residualForcesMovingAverageHistorySample.pop();
// "Time/(#nodes*#iterations)",
// timePerNodePerIterationHistory);
}
if ((mSettings.shouldCreatePlots || mSettings.isDebugMode) && mCurrentSimulationStep != 0) { pMesh->residualForcesMovingAverageDerivativeNorm
history.stepPulse(*pMesh); = std::abs((residualForcesMovingAverageHistorySample.back()
} - residualForcesMovingAverageHistorySample.front()))
/ (movingAverageSampleSize);
residualForcesMovingAverageDerivativeMax
= std::max(pMesh->residualForcesMovingAverageDerivativeNorm,
residualForcesMovingAverageDerivativeMax);
pMesh->residualForcesMovingAverageDerivativeNorm
/= residualForcesMovingAverageDerivativeMax;
// std::cout << "Normalized derivative:"
// << residualForcesMovingAverageDerivativeNorm
// << std::endl;
}
if (std::isnan(pMesh->currentTotalKineticEnergy)) { // pMesh->previousTotalPotentialEnergykN =
if (mSettings.beVerbose) { // pMesh->currentTotalPotentialEnergykN;
std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl; // pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000;
} // timePerNodePerIterationHistor.push_back(timePerNodePerIteration);
break; if (mSettings.beVerbose && mSettings.isDebugMode
} && mCurrentSimulationStep % mSettings.debugModeStep == 0) {
printCurrentState();
std::cout << "Percentage of target:"
<< 100 * mSettings.totalResidualForcesNormThreshold
/ pMesh->totalResidualForcesNorm
<< "%" << std::endl;
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver",
// movingAverages);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver deriv",
// movingAveragesDerivatives);
// draw();
// SimulationResulnodalForcedDisplacementstsReporter::createPlot("Number of Steps",
// "Time/(#nodes*#iterations)",
// timePerNodePerIterationHistory);
}
if ((mSettings.shouldCreatePlots || mSettings.isDebugMode) && mCurrentSimulationStep != 0) {
history.stepPulse(*pMesh);
percentageOfConvergence.push_back(100 * mSettings.totalResidualForcesNormThreshold
/ pMesh->totalResidualForcesNorm);
}
if (std::isnan(pMesh->currentTotalKineticEnergy)) {
if (mSettings.beVerbose) {
std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl;
}
break;
}
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if ((mSettings.shouldDraw && mCurrentSimulationStep % mSettings.drawingStep == 0) if (mSettings.shouldDraw && mSettings.isDebugMode
|| (mSettings.isDebugMode && mCurrentSimulationStep % mSettings.debugModeStep == 0)) /* && && mCurrentSimulationStep % mSettings.drawingStep == 0) /* &&
currentSimulationStep > maxDRMIterations*/ currentSimulationStep > maxDRMIterations*/
{ {
// std::string saveTo = std::filesystem::current_path() // std::string saveTo = std::filesystem::current_path()
// .append("Debugging_files") // .append("Debugging_files")
// .append("Screenshots") // .append("Screenshots")
// .string(); // .string();
draw(); draw();
// yValues = history.kineticEnergy; // yValues = history.kineticEnergy;
// plotHandle = matplot::scatter(xPlot, yValues); // plotHandle = matplot::scatter(xPlot, yValues);
// label = "Log of Kinetic energy"; // label = "Log of Kinetic energy";
// plotHandle->legend_string(label); // plotHandle->legend_string(label);
// shouldUseKineticEnergyThreshold = true; // shouldUseKineticEnergyThreshold = true;
} }
#endif #endif
if (mSettings.isDebugMode && mCurrentSimulationStep % mSettings.debugModeStep == 0) { if (mSettings.shouldCreatePlots && mSettings.isDebugMode
&& mCurrentSimulationStep % mSettings.debugModeStep == 0) {
// SimulationResultsReporter::createPlot("Number of Steps", // SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver deriv", // "Residual Forces mov aver deriv",
// movingAveragesDerivatives_norm); // movingAveragesDerivatives_norm);
SimulationResultsReporter::createPlot("Number of Steps",
"Residual Forces mov aver",
history.residualForcesMovingAverage);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces",
// history.residualForces);
// SimulationResultsReporter::createPlot("Number of Steps", // SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Kinetic energy", // "Residual Forces mov aver",
// history.kineticEnergy); // history.residualForcesMovingAverage);
SimulationResultsReporter::createPlot("Number of Steps",
"Log of Residual Forces",
history.logResidualForces,
{},
history.redMarks);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Kinetic energy",
// history.kineticEnergy);
// SimulationResultsReporter reporter;
// reporter.reportHistory(history, "Results");
// SimulationResultsReporter::createPlot("Number of Steps",
// "Percentage of convergence",
// percentageOfConvergence,
// {},
// history.redMarks);
} }
// t = t + Dt; // t = t + Dt;
// std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy // std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy
@ -1990,16 +2206,19 @@ mesh->currentTotalPotentialEnergykN*/
= mSettings.shouldUseTranslationalKineticEnergyThreshold = mSettings.shouldUseTranslationalKineticEnergyThreshold
&& pMesh->currentTotalTranslationalKineticEnergy && pMesh->currentTotalTranslationalKineticEnergy
< mSettings.totalTranslationalKineticEnergyThreshold < mSettings.totalTranslationalKineticEnergyThreshold
&& mCurrentSimulationStep > 20 && numOfDampings > 2; && mCurrentSimulationStep > 20 && numOfDampings > 0;
const bool fullfillsResidualForcesNormTerminationCriterion const bool fullfillsResidualForcesNormTerminationCriterion
= pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold; = pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold;
const bool fullfillsMovingAverageNormTerminationCriterion const bool fullfillsMovingAverageNormTerminationCriterion
= pMesh->residualForcesMovingAverage = pMesh->residualForcesMovingAverage
< mSettings.residualForcesMovingAverageNormThreshold; < mSettings.residualForcesMovingAverageNormThreshold;
/*pMesh->residualForcesMovingAverage < totalResidualForcesNormThreshold*/ /*pMesh->residualForcesMovingAverage < totalResidualForcesNormThreshold*/
const bool fullfillsMovingAverageDerivativeNormTerminationCriterion
= pMesh->residualForcesMovingAverageDerivativeNorm < 1e-4;
const bool shouldTerminate const bool shouldTerminate
= fullfillsKineticEnergyTerminationCriterion = fullfillsKineticEnergyTerminationCriterion
// || fullfillsMovingAverageNormTerminationCriterion // || fullfillsMovingAverageNormTerminationCriterion
// || fullfillsMovingAverageDerivativeNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion; || fullfillsResidualForcesNormTerminationCriterion;
if (shouldTerminate) { if (shouldTerminate) {
if (mSettings.beVerbose) { if (mSettings.beVerbose) {
@ -2021,34 +2240,59 @@ mesh->currentTotalPotentialEnergykN*/
break; break;
// } // }
} }
// printCurrentState();
// for (VertexType &v : mesh->vert) {
// Node &node = mesh->nodes[v];
// node.displacements = node.displacements - node.velocity * Dt;
// }
// applyDisplacements(constrainedVertices);
// if (!pJob->nodalForcedDisplacements.empty()) {
// applyForcedDisplacements(
// pJob->nodalForcedDisplacements); // if (mSettings.beVerbose) {
// } // printCurrentState();
// updateElementalLengths(); // }
// 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();
// const double triggerPercentage = 0.001;
// const double xi_min = 0.85;
// const double xi_init = 0.9969;
// if (mSettings.totalResidualForcesNormThreshold / pMesh->totalResidualForcesNorm
// > triggerPercentage) {
// mSettings.xi = ((xi_min - xi_init)
// * (mSettings.totalResidualForcesNormThreshold
// / pMesh->totalResidualForcesNorm)
// + xi_init - triggerPercentage * xi_min)
// / (1 - triggerPercentage);
// }
Dt *= mSettings.xi;
resetVelocities(); resetVelocities();
Dt = Dt * mSettings.xi;
++numOfDampings; ++numOfDampings;
// std::cout << "Number of dampings:" << numOfDampings << endl; if (mSettings.shouldCreatePlots) {
history.redMarks.push_back(mCurrentSimulationStep);
}
// std::cout << "Number of dampings:" << numOfDampings << endl;
// draw();
} }
} }
SimulationResults results = computeResults(pJob); SimulationResults results = computeResults(pJob);
if (mCurrentSimulationStep == settings.maxDRMIterations && if (mCurrentSimulationStep == settings.maxDRMIterations && mCurrentSimulationStep != 0) {
mCurrentSimulationStep != 0) { std::cout << "Did not reach equilibrium before reaching the maximum number "
std::cout << "Did not reach equilibrium before reaching the maximum number " "of DRM steps ("
"of DRM steps (" << settings.maxDRMIterations << "). Breaking simulation" << std::endl;
<< settings.maxDRMIterations << "). Breaking simulation" results.converged = false;
<< std::endl;
results.converged = false;
} else if (std::isnan(pMesh->currentTotalKineticEnergy)) { } else if (std::isnan(pMesh->currentTotalKineticEnergy)) {
std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl; std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl;
results.converged = false; results.converged = false;
@ -2064,7 +2308,6 @@ mesh->currentTotalPotentialEnergykN*/
<< "s" << std::endl; << "s" << std::endl;
std::cout << "Number of dampings:" << numOfDampings << endl; std::cout << "Number of dampings:" << numOfDampings << endl;
std::cout << "Number of steps:" << mCurrentSimulationStep << endl; std::cout << "Number of steps:" << mCurrentSimulationStep << endl;
results.converged = true;
} }
// mesh.printVertexCoordinates(mesh.VN() / 2); // mesh.printVertexCoordinates(mesh.VN() / 2);
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
@ -2073,46 +2316,93 @@ mesh->currentTotalPotentialEnergykN*/
} }
#endif #endif
if (results.converged) { if (results.converged) {
results.debug_drmDisplacements = results.displacements;
results.internalPotentialEnergy = computeTotalInternalPotentialEnergy(); results.internalPotentialEnergy = computeTotalInternalPotentialEnergy();
results.rotationalDisplacementQuaternion.resize(pMesh->VN()); results.rotationalDisplacementQuaternion.resize(pMesh->VN());
results.debug_q_f1.resize(pMesh->VN());
results.debug_q_normal.resize(pMesh->VN());
results.debug_q_nr.resize(pMesh->VN());
for (int vi = 0; vi < pMesh->VN(); vi++) { for (int vi = 0; vi < pMesh->VN(); vi++) {
const Node &node = pMesh->nodes[vi]; const Node &node = pMesh->nodes[vi];
const Eigen::Vector3d initialNormal = node.initialNormal.ToEigenVector<Eigen::Vector3d>(); const Eigen::Vector3d nInitial_eigen = node.initialNormal.ToEigenVector<Eigen::Vector3d>();
const Eigen::Vector3d deformedNormal const Eigen::Vector3d nDeformed_eigen
= pMesh->vert[vi].cN().ToEigenVector<Eigen::Vector3d>(); = pMesh->vert[vi].cN().ToEigenVector<Eigen::Vector3d>();
Eigen::Quaternion<double> q; Eigen::Quaternion<double> q_normal;
q.setFromTwoVectors(initialNormal, deformedNormal); q_normal.setFromTwoVectors(nInitial_eigen, nDeformed_eigen);
Eigen::Quaternion<double> q_nr; Eigen::Quaternion<double> q_nr_nDeformed;
q_nr = Eigen::AngleAxis<double>(pMesh->nodes[vi].nR, deformedNormal); q_nr_nDeformed = Eigen::AngleAxis<double>(pMesh->nodes[vi].nR, nDeformed_eigen);
Eigen::Quaternion<double> q_nr_nInit;
q_nr_nInit = Eigen::AngleAxis<double>(pMesh->nodes[vi].nR, nInitial_eigen);
const auto w = q_nr_nDeformed.w();
const auto theta = 2 * acos(q_nr_nDeformed.w());
const auto nr = pMesh->nodes[vi].nR;
Eigen::Vector3d deformedNormal_debug(q_nr_nDeformed.x() * sin(theta / 2),
q_nr_nDeformed.y() * sin(theta / 2),
q_nr_nDeformed.z() * sin(theta / 2));
deformedNormal_debug.normalize();
// const double nr = nr_0To2pi <= M_PI ? nr_0To2pi : (nr_0To2pi - 2 * M_PI);
const double nr_debug = deformedNormal_debug.dot(nDeformed_eigen) > 0 ? theta : -theta;
assert(pMesh->nodes[vi].nR - nr_debug < 1e-6);
VectorType referenceT1_deformed = pMesh->elements[node.referenceElement].frame.t1; VectorType referenceT1_deformed = pMesh->elements[node.referenceElement].frame.t1;
const VectorType &n_deformed = pMesh->vert[vi].cN(); const VectorType &nDeformed = pMesh->vert[vi].cN();
const VectorType referenceF1_deformed = (referenceT1_deformed const VectorType referenceF1_deformed = (referenceT1_deformed
- (n_deformed - (node.initialNormal
* (referenceT1_deformed * n_deformed))) * (referenceT1_deformed * node.initialNormal)))
.Normalize(); .Normalize();
const VectorType referenceT1_initial const VectorType referenceT1_initial
= computeT1Vector(pMesh->nodes[node.referenceElement->cV(0)].initialLocation, = computeT1Vector(pMesh->nodes[node.referenceElement->cV(0)].initialLocation,
pMesh->nodes[node.referenceElement->cV(1)].initialLocation); pMesh->nodes[node.referenceElement->cV(1)].initialLocation);
const VectorType &n_initial = node.initialNormal; // const VectorType &n_initial = node.initialNormal;
const VectorType referenceF1_initial const VectorType referenceF1_initial = (referenceT1_initial
= (referenceT1_initial - (n_deformed * (referenceT1_initial * n_deformed))).Normalize(); - (node.initialNormal
Eigen::Quaternion<double> q_f1; //nr is with respect to f1 * (referenceT1_initial * node.initialNormal)))
q_f1.setFromTwoVectors(referenceF1_initial.ToEigenVector<Eigen::Vector3d>(), .Normalize();
referenceF1_deformed.ToEigenVector<Eigen::Vector3d>()); Eigen::Quaternion<double> q_f1_nInit; //nr is with respect to f1
Eigen::Quaternion<double> q_t1; q_f1_nInit.setFromTwoVectors(referenceF1_initial.ToEigenVector<Eigen::Vector3d>(),
q_t1.setFromTwoVectors(referenceT1_initial.ToEigenVector<Eigen::Vector3d>(), referenceF1_deformed.ToEigenVector<Eigen::Vector3d>());
referenceT1_deformed.ToEigenVector<Eigen::Vector3d>());
results.rotationalDisplacementQuaternion[vi] = q * q_f1 * q_nr; Eigen::Quaternion<double> q_f1_nDeformed; //nr is with respect to f1
// const VectorType &n_initial = node.initialNormal;
const VectorType referenceF1_initial_def
= (referenceT1_initial - (nDeformed * (referenceT1_initial * nDeformed))).Normalize();
const VectorType referenceF1_deformed_def
= (referenceT1_deformed - (nDeformed * (referenceT1_deformed * nDeformed))).Normalize();
q_f1_nDeformed.setFromTwoVectors(referenceF1_initial_def.ToEigenVector<Eigen::Vector3d>(),
referenceF1_deformed_def.ToEigenVector<Eigen::Vector3d>());
const auto debug_qf1_nDef = (q_f1_nDeformed * q_nr_nDeformed) * nDeformed_eigen;
const auto debug_qf1_nInit = (q_f1_nInit * q_nr_nInit) * nInitial_eigen;
const auto debug_deformedNormal_f1Init = (q_normal * (q_f1_nInit * q_nr_nInit))
* nInitial_eigen;
const auto debug_deformedNormal_f1Def = ((q_nr_nDeformed * q_f1_nDeformed) * q_normal)
* nInitial_eigen;
// Eigen::Quaternion<double> q_t1;
// q_t1.setFromTwoVectors(referenceT1_initial.ToEigenVector<Eigen::Vector3d>(),
// referenceT1_deformed.ToEigenVector<Eigen::Vector3d>());
results.debug_q_f1[vi] = q_f1_nInit;
results.debug_q_normal[vi] = q_normal;
results.debug_q_nr[vi] = q_nr_nInit;
results.rotationalDisplacementQuaternion[vi]
= (q_normal * (q_f1_nInit * q_nr_nInit)); //q_f1_nDeformed * q_nr_nDeformed * q_normal;
//Update the displacement vector to contain the euler angles //Update the displacement vector to contain the euler angles
const Eigen::Vector3d eulerAngles const Eigen::Vector3d eulerAngles
= results.rotationalDisplacementQuaternion[vi].toRotationMatrix().eulerAngles(0, 1, 2); = results.rotationalDisplacementQuaternion[vi].toRotationMatrix().eulerAngles(0, 1, 2);
results.displacements[vi][3] = eulerAngles[0]; results.displacements[vi][3] = eulerAngles[0];
results.displacements[vi][4] = eulerAngles[1]; results.displacements[vi][4] = eulerAngles[1];
results.displacements[vi][5] = eulerAngles[2]; results.displacements[vi][5] = eulerAngles[2];
// Eigen::Quaterniond q_test = Eigen::AngleAxisd(eulerAngles[0], Eigen::Vector3d::UnitX())
// * Eigen::AngleAxisd(eulerAngles[1], Eigen::Vector3d::UnitY())
// * Eigen::AngleAxisd(eulerAngles[2], Eigen::Vector3d::UnitZ());
// const double dot_test = results.rotationalDisplacementQuaternion[vi].dot(q_test);
// assert(dot_test > 1 - 1e5);
int i = 0;
i++;
} }
} }

122
drmsimulationmodel.hpp
View File

@ -26,7 +26,7 @@ public:
struct Settings struct Settings
{ {
bool isDebugMode{false}; bool isDebugMode{false};
int debugModeStep{10000}; int debugModeStep{100000};
bool shouldDraw{false}; bool shouldDraw{false};
bool beVerbose{false}; bool beVerbose{false};
bool shouldCreatePlots{false}; bool shouldCreatePlots{false};
@ -40,79 +40,78 @@ public:
int maxDRMIterations{0}; int maxDRMIterations{0};
bool shouldUseTranslationalKineticEnergyThreshold{false}; bool shouldUseTranslationalKineticEnergyThreshold{false};
int gradualForcedDisplacementSteps{100}; int gradualForcedDisplacementSteps{100};
double gamma{0.8};
std::optional<double> displacementCap;
double totalResidualForcesNormThreshold{1e-3};
Settings() {} Settings() {}
}; };
double totalResidualForcesNormThreshold{1};
private: private:
Settings mSettings; Settings mSettings;
double Dt{mSettings.Dtini}; double Dt{mSettings.Dtini};
bool checkedForMaximumMoment{false}; bool checkedForMaximumMoment{false};
double externalMomentsNorm{0}; bool shouldTemporarilyDampForces{false};
size_t mCurrentSimulationStep{0}; bool shouldTemporarilyAmplifyForces{true};
matplot::line_handle plotHandle; double externalMomentsNorm{0};
std::vector<double> plotYValues; size_t mCurrentSimulationStep{0};
size_t numOfDampings{0}; matplot::line_handle plotHandle;
std::vector<double> plotYValues;
size_t numOfDampings{0};
const std::string meshPolyscopeLabel{"Simulation mesh"}; const std::string meshPolyscopeLabel{"Simulation mesh"};
std::unique_ptr<SimulationMesh> pMesh; std::unique_ptr<SimulationMesh> pMesh;
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> std::unordered_map<VertexIndex, std::unordered_set<DoFType>> constrainedVertices;
constrainedVertices; 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; std::unordered_set<VertexIndex> rigidSupports;
std::unordered_set<VertexIndex> rigidSupports;
void reset(); void reset();
void updateNodalInternalForces( void updateNodalInternalForces(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
&fixedVertices); void updateNodalExternalForces(
void updateNodalExternalForces( const std::unordered_map<VertexIndex, Vector6d> &nodalForces,
const std::unordered_map<VertexIndex, Vector6d> &nodalForces, const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> void updateResidualForces();
&fixedVertices); void updateRotationalDisplacements();
void updateResidualForces(); void updateElementalLengths();
void updateRotationalDisplacements();
void updateElementalLengths();
void updateNodalMasses(); void updateNodalMasses();
void updateNodalAccelerations(); void updateNodalAccelerations();
void updateNodalVelocities(); void updateNodalVelocities();
void updateNodalDisplacements(); void updateNodalDisplacements();
void applyForcedDisplacements( void applyForcedDisplacements(
const std::unordered_map<VertexIndex, Eigen::Vector3d> &nodalForcedDisplacements); const std::unordered_map<VertexIndex, Eigen::Vector3d> &nodalForcedDisplacements);
Vector6d computeElementTorsionalForce( Vector6d computeElementTorsionalForce(const Element &element,
const Element &element, const Vector6d &displacementDifference, const Vector6d &displacementDifference,
const std::unordered_set<DoFType> &constrainedDof); const std::unordered_set<DoFType> &constrainedDof);
// BeamFormFinder::Vector6d computeElementFirstBendingForce( // BeamFormFinder::Vector6d computeElementFirstBendingForce(
// const Element &element, const Vector6d &displacementDifference, // const Element &element, const Vector6d &displacementDifference,
// const std::unordered_set<gsl::index> &constrainedDof); // const std::unordered_set<gsl::index> &constrainedDof);
// BeamFormFinder::Vector6d computeElementSecondBendingForce( // BeamFormFinder::Vector6d computeElementSecondBendingForce(
// const Element &element, const Vector6d &displacementDifference, // const Element &element, const Vector6d &displacementDifference,
// const std::unordered_set<gsl::index> &constrainedDof); // const std::unordered_set<gsl::index> &constrainedDof);
void updateKineticEnergy(); void updateKineticEnergy();
void resetVelocities(); void resetVelocities();
SimulationResults computeResults(const std::shared_ptr<SimulationJob> &pJob); SimulationResults computeResults(const std::shared_ptr<SimulationJob> &pJob);
void updateNodePosition( void updateNodePosition(
VertexType &v, VertexType &v,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
&fixedVertices);
void applyDisplacements( void applyDisplacements(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
&fixedVertices);
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
void draw(const string &screenshotsFolder= {}); void draw(const string &screenshotsFolder= {});
@ -206,11 +205,16 @@ public:
double computeTotalInternalPotentialEnergy(); double computeTotalInternalPotentialEnergy();
void applySolutionGuess(const SimulationResults &solutionGuess,
const std::shared_ptr<SimulationJob> &pJob);
void updateNodeNr(VertexType &v);
public: public:
DRMSimulationModel(); DRMSimulationModel();
SimulationResults SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
executeSimulation(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings = Settings(),
const Settings &settings = Settings()); const SimulationResults &solutionGuess = SimulationResults());
static void runUnitTests(); static void runUnitTests();
}; };

1
edgemesh.cpp
View File

@ -220,6 +220,7 @@ bool VCGEdgeMesh::loadUsingNanoply(const std::string &plyFilename) {
unsigned int mask = 0; unsigned int mask = 0;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD; mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL; mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOLOR;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX; mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::LoadModel(plyFilename.c_str(), if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::LoadModel(plyFilename.c_str(),
*this, mask) != 0) { *this, mask) != 0) {

7
polymesh.hpp
View File

@ -110,7 +110,7 @@ public:
return faces; return faces;
} }
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
std::shared_ptr<polyscope::SurfaceMesh> registerForDrawing( polyscope::SurfaceMesh *registerForDrawing(
const std::optional<std::array<double, 3>> &desiredColor = std::nullopt, const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
const bool &shouldEnable = true) const bool &shouldEnable = true)
{ {
@ -118,8 +118,9 @@ public:
auto faces = getFaces(); auto faces = getFaces();
PolyscopeInterface::init(); PolyscopeInterface::init();
std::shared_ptr<polyscope::SurfaceMesh> polyscopeHandle_mesh( polyscope::SurfaceMesh *polyscopeHandle_mesh = polyscope::registerSurfaceMesh(label,
polyscope::registerSurfaceMesh(label, vertices, faces)); vertices,
faces);
const double drawingRadius = 0.002; const double drawingRadius = 0.002;
polyscopeHandle_mesh->setEnabled(shouldEnable); polyscopeHandle_mesh->setEnabled(shouldEnable);

64
reducedmodeloptimizer_structs.hpp
View File

@ -362,7 +362,7 @@ struct Colors
#endif #endif
} }
bool load(const string &loadFromPath) bool load(const string &loadFromPath, const bool &shouldLoadDebugFiles = false)
{ {
assert(std::filesystem::is_directory(loadFromPath)); assert(std::filesystem::is_directory(loadFromPath));
std::filesystem::path jsonFilepath( std::filesystem::path jsonFilepath(
@ -414,37 +414,39 @@ struct Colors
baseTriangle.P2(0) = CoordType(baseTriangleVertices[6], baseTriangle.P2(0) = CoordType(baseTriangleVertices[6],
baseTriangleVertices[7], baseTriangleVertices[7],
baseTriangleVertices[8]); baseTriangleVertices[8]);
if (shouldLoadDebugFiles) {
const std::filesystem::path simulationJobsFolderPath( const std::filesystem::path simulationJobsFolderPath(
std::filesystem::path(loadFromPath).append("SimulationJobs")); std::filesystem::path(loadFromPath).append("SimulationJobs"));
for (const auto &directoryEntry : for (const auto &directoryEntry :
filesystem::directory_iterator(simulationJobsFolderPath)) { filesystem::directory_iterator(simulationJobsFolderPath)) {
const auto simulationScenarioPath = directoryEntry.path(); const auto simulationScenarioPath = directoryEntry.path();
if (!std::filesystem::is_directory(simulationScenarioPath)) { if (!std::filesystem::is_directory(simulationScenarioPath)) {
continue; continue;
}
// Load full pattern files
for (const auto &fileEntry : filesystem::directory_iterator(
std::filesystem::path(simulationScenarioPath).append("Full"))) {
const auto filepath = fileEntry.path();
if (filepath.extension() == ".json") {
SimulationJob job;
job.load(filepath.string());
fullPatternSimulationJobs.push_back(std::make_shared<SimulationJob>(job));
} }
} // Load full pattern files
for (const auto &fileEntry : filesystem::directory_iterator(
std::filesystem::path(simulationScenarioPath).append("Full"))) {
const auto filepath = fileEntry.path();
if (filepath.extension() == ".json") {
SimulationJob job;
// Load reduced pattern files and apply optimal parameters job.load(filepath.string());
for (const auto &fileEntry : filesystem::directory_iterator( fullPatternSimulationJobs.push_back(std::make_shared<SimulationJob>(job));
std::filesystem::path(simulationScenarioPath).append("Reduced"))) { }
const auto filepath = fileEntry.path(); }
if (filepath.extension() == ".json") {
SimulationJob job; // Load reduced pattern files and apply optimal parameters
job.load(filepath.string()); for (const auto &fileEntry : filesystem::directory_iterator(
applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh); std::filesystem::path(simulationScenarioPath).append("Reduced"))) {
applyOptimizationResults_elements(*this, job.pMesh); const auto filepath = fileEntry.path();
reducedPatternSimulationJobs.push_back(std::make_shared<SimulationJob>(job)); if (filepath.extension() == ".json") {
SimulationJob job;
job.load(filepath.string());
applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh);
applyOptimizationResults_elements(*this, job.pMesh);
reducedPatternSimulationJobs.push_back(
std::make_shared<SimulationJob>(job));
}
} }
} }
} }
@ -575,7 +577,7 @@ struct Colors
SimulationResults reducedModelLinearResults = linearSimulator.executeSimulation( SimulationResults reducedModelLinearResults = linearSimulator.executeSimulation(
pReducedPatternSimulationJob); pReducedPatternSimulationJob);
reducedModelLinearResults.setLabelPrefix("linear"); reducedModelLinearResults.setLabelPrefix("linear");
reducedModelLinearResults.registerForDrawing(Colors::reducedDeformed, true, true); reducedModelLinearResults.registerForDrawing(Colors::reducedDeformed, false);
polyscope::options::programName = fullPatternSimulationJobs[0]->pMesh->getLabel(); polyscope::options::programName = fullPatternSimulationJobs[0]->pMesh->getLabel();
polyscope::view::resetCameraToHomeView(); polyscope::view::resetCameraToHomeView();
polyscope::show(); polyscope::show();

48
simulation_structs.hpp
View File

@ -35,7 +35,7 @@ struct SimulationHistory {
size_t numberOfSteps{0}; size_t numberOfSteps{0};
std::string label; std::string label;
std::vector<double> residualForcesLog; std::vector<double> logResidualForces;
std::vector<double> kineticEnergy; std::vector<double> kineticEnergy;
std::vector<double> potentialEnergies; std::vector<double> potentialEnergies;
std::vector<size_t> redMarks; std::vector<size_t> redMarks;
@ -48,19 +48,20 @@ struct SimulationHistory {
void markGreen(const size_t &stepNumber) { greenMarks.push_back(stepNumber); } void markGreen(const size_t &stepNumber) { greenMarks.push_back(stepNumber); }
void stepPulse(const SimulationMesh &mesh) { void stepPulse(const SimulationMesh &mesh)
kineticEnergy.push_back(log(mesh.currentTotalKineticEnergy)); {
// potentialEnergy.push_back(mesh.totalPotentialEnergykN); kineticEnergy.push_back(std::log10(mesh.currentTotalKineticEnergy));
residualForcesLog.push_back(std::log(mesh.totalResidualForcesNorm)); // potentialEnergy.push_back(mesh.totalPotentialEnergykN);
residualForcesMovingAverage.push_back(mesh.residualForcesMovingAverage); logResidualForces.push_back(std::log10(mesh.totalResidualForcesNorm));
sumOfNormalizedDisplacementNorms.push_back(mesh.sumOfNormalizedDisplacementNorms); residualForcesMovingAverage.push_back(mesh.residualForcesMovingAverage);
// residualForcesMovingAverageDerivativesLog.push_back( sumOfNormalizedDisplacementNorms.push_back(mesh.sumOfNormalizedDisplacementNorms);
// std::log(mesh.residualForcesMovingAverageDerivativeNorm)); // residualForcesMovingAverageDerivativesLog.push_back(
// std::log(mesh.residualForcesMovingAverageDerivativeNorm));
} }
void clear() void clear()
{ {
residualForcesLog.clear(); logResidualForces.clear();
kineticEnergy.clear(); kineticEnergy.clear();
potentialEnergies.clear(); potentialEnergies.clear();
residualForcesMovingAverage.clear(); residualForcesMovingAverage.clear();
@ -119,6 +120,11 @@ public:
SimulationJob() {} SimulationJob() {}
SimulationJob(const std::string &jsonFilename) { load(jsonFilename); } SimulationJob(const std::string &jsonFilename) { load(jsonFilename); }
bool isEmpty()
{
return constrainedVertices.empty() && nodalExternalForces.empty()
&& nodalForcedDisplacements.empty() && pMesh == nullptr;
}
void remap(const std::unordered_map<size_t, size_t> &thisToOtherViMap, void remap(const std::unordered_map<size_t, size_t> &thisToOtherViMap,
SimulationJob &simulationJobMapped) const SimulationJob &simulationJobMapped) const
{ {
@ -418,14 +424,18 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
->setEnabled(shouldEnable); ->setEnabled(shouldEnable);
} }
// per node external moments // per node external moments
bool hasExternalMoments = false;
std::vector<std::array<double, 3>> externalMoments(pMesh->VN()); std::vector<std::array<double, 3>> externalMoments(pMesh->VN());
for (const auto &forcePair : nodalExternalForces) { for (const auto &forcePair : nodalExternalForces) {
auto index = forcePair.first; auto index = forcePair.first;
auto load = forcePair.second; const Vector6d &load = forcePair.second;
if (load.getRotation().norm() != 0) {
hasExternalMoments = true;
}
externalMoments[index] = {load[3], load[4], load[5]}; externalMoments[index] = {load[3], load[4], load[5]};
} }
if (!externalMoments.empty()) { if (hasExternalMoments) {
polyscope::getCurveNetwork(meshLabel) polyscope::getCurveNetwork(meshLabel)
->addNodeVectorQuantity("External moment_" + label, externalMoments) ->addNodeVectorQuantity("External moment_" + label, externalMoments)
->setEnabled(shouldEnable); ->setEnabled(shouldEnable);
@ -436,8 +446,12 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
if (polyscope::getCurveNetwork(meshLabel) == nullptr) { if (polyscope::getCurveNetwork(meshLabel) == nullptr) {
return; return;
} }
polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force_" + label); if (!nodalExternalForces.empty()) {
polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label); polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force_" + label);
}
if (!constrainedVertices.empty()) {
polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label);
}
} }
#endif // POLYSCOPE_DEFINED #endif // POLYSCOPE_DEFINED
}; };
@ -446,9 +460,13 @@ struct SimulationResults
/*TODO: remove rotationalDisplacementQuaternion since the last three components of the displacments /*TODO: remove rotationalDisplacementQuaternion since the last three components of the displacments
* vector contains the same info using euler angles * vector contains the same info using euler angles
*/ */
bool converged{true}; bool converged{false};
std::shared_ptr<SimulationJob> job; std::shared_ptr<SimulationJob> job;
SimulationHistory history; SimulationHistory history;
std::vector<Vector6d> debug_drmDisplacements;
std::vector<Eigen::Quaternion<double>> debug_q_f1; //per vertex
std::vector<Eigen::Quaternion<double>> debug_q_normal; //per vertex
std::vector<Eigen::Quaternion<double>> debug_q_nr; //per vertex
std::vector<Vector6d> displacements; std::vector<Vector6d> displacements;
std::vector<Eigen::Quaternion<double>> rotationalDisplacementQuaternion; //per vertex std::vector<Eigen::Quaternion<double>> rotationalDisplacementQuaternion; //per vertex
double internalPotentialEnergy{0}; double internalPotentialEnergy{0};

85
simulationhistoryplotter.hpp
View File

@ -38,13 +38,13 @@ struct SimulationResultsReporter {
file << "\n"; file << "\n";
} }
if (!history.residualForcesLog.empty()) { if (!history.logResidualForces.empty()) {
file << "Residual forces" file << "Residual forces"
<< "\n"; << "\n";
for (size_t step = 0; step < numberOfSteps; step++) { for (size_t step = 0; step < numberOfSteps; step++) {
file << history.residualForcesLog[step] << "\n"; file << history.logResidualForces[step] << "\n";
} }
file << "\n"; file << "\n";
} }
if (!history.potentialEnergies.empty()) { if (!history.potentialEnergies.empty()) {
@ -58,6 +58,15 @@ struct SimulationResultsReporter {
file.close(); file.close();
} }
void reportHistory(const SimulationHistory &history,
const std::string &reportFolderPath,
const std::string &graphSuffix = std::string())
{
const auto simulationResultPath = std::filesystem::path(reportFolderPath).append(history.label);
std::filesystem::create_directory(simulationResultPath);
createPlots(history, simulationResultPath, graphSuffix);
}
void reportResults(const std::vector<SimulationResults> &results, void reportResults(const std::vector<SimulationResults> &results,
const std::string &reportFolderPath, const std::string &reportFolderPath,
const std::string &graphSuffix = std::string()) { const std::string &graphSuffix = std::string()) {
@ -78,45 +87,55 @@ struct SimulationResultsReporter {
writeStatistics(simulationResult, simulationResultPath.string()); writeStatistics(simulationResult, simulationResultPath.string());
} }
} }
static void createPlot(const std::string &xLabel, const std::string &yLabel, static void createPlot(const std::string &xLabel,
const std::vector<double> &XvaluesToPlot, const std::string &yLabel,
const std::vector<double> &YvaluesToPlot, const std::vector<double> &x,
const std::string &saveTo = {}) { const std::vector<double> &y,
// matplot::figure(true); const std::vector<double> &markerSizes,
matplot::xlabel(xLabel); const std::vector<double> &c,
matplot::ylabel(yLabel); const std::string &saveTo = {})
matplot::grid(matplot::on); {
matplot::scatter(XvaluesToPlot, YvaluesToPlot); // matplot::figure(true);
if (!saveTo.empty()) { matplot::xlabel(xLabel);
matplot::save(saveTo); matplot::ylabel(yLabel);
} matplot::grid(matplot::on);
matplot::scatter(x, y, markerSizes, c)->marker_face(true);
if (!saveTo.empty()) {
matplot::save(saveTo);
}
} }
static void createPlot(const std::string &xLabel, static void createPlot(const std::string &xLabel,
const std::string &yLabel, const std::string &yLabel,
const std::vector<double> &YvaluesToPlot, const std::vector<double> &YvaluesToPlot,
const std::string &saveTo = {}) const std::string &saveTo = {},
const std::vector<size_t> &markPoints = {})
{ {
if (YvaluesToPlot.size() < 2) { if (YvaluesToPlot.size() < 2) {
return; return;
} }
auto x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size()); auto x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size());
createPlot(xLabel, yLabel, x, YvaluesToPlot, saveTo); std::vector<double> colors(x.size(), 0.5);
// matplot::figure(true); std::vector<double> markerSizes(x.size(), 5);
// matplot::hold(matplot::on); if (!markPoints.empty()) {
for (const auto pointIndex : markPoints) {
colors[pointIndex] = 0.9;
markerSizes[pointIndex] = 14;
}
}
createPlot(xLabel, yLabel, x, YvaluesToPlot, markerSizes, colors, saveTo);
// matplot::figure(true);
// matplot::hold(matplot::on);
// ->marker_indices(history.redMarks) // ->marker_indices(history.redMarks)
// ->marker_indices(truncatedRedMarks) // ->marker_indices(truncatedRedMarks)
// .marker_color("r") // .marker_color("r")
// ->marker_size(1) // ->marker_size(1)
; // ;
// auto greenMarksY = matplot::transform( // auto greenMarksY = matplot::transform(history.greenMarks,
// history.greenMarks, [&](auto x) { return history.kineticEnergy[x]; // [&](auto x) { return YvaluesToPlot[x]; });
// }); // matplot::scatter(history.greenMarks, greenMarksY)->color("green").marker_size(10);
// matplot::scatter(history.greenMarks, greenMarksY) // matplot::hold(matplot::off);
// ->color("green")
// .marker_size(10);
// matplot::hold(matplot::off);
} }
void createPlots(const SimulationHistory &history, void createPlots(const SimulationHistory &history,
@ -136,10 +155,10 @@ struct SimulationResultsReporter {
.string()); .string());
} }
if (!history.residualForcesLog.empty()) { if (!history.logResidualForces.empty()) {
createPlot("Number of Iterations", createPlot("Number of Iterations",
"Residual Forces norm log", "Residual Forces norm log",
history.residualForcesLog, history.logResidualForces,
std::filesystem::path(graphsFolder) std::filesystem::path(graphsFolder)
.append("ResidualForcesLog_" + graphSuffix + ".png") .append("ResidualForcesLog_" + graphSuffix + ".png")
.string()); .string());

68
simulationmesh.cpp
View File

@ -40,20 +40,46 @@ SimulationMesh::SimulationMesh(PatternGeometry &pattern) {
initializeElements(); initializeElements();
} }
SimulationMesh::SimulationMesh(SimulationMesh &mesh) { SimulationMesh::SimulationMesh(SimulationMesh &m)
vcg::tri::MeshAssert<SimulationMesh>::VertexNormalNormalized(mesh); {
VCGEdgeMesh::copy(mesh); vcg::tri::MeshAssert<SimulationMesh>::VertexNormalNormalized(m);
VCGEdgeMesh::copy(m);
elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>( elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(*this,
*this, std::string("Elements")); std::string(
nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>( "Elements"));
*this, std::string("Nodes")); nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(*this,
initializeNodes(); std::string("Nodes"));
initializeNodes();
for (size_t ei = 0; ei < EN(); ei++) { for (size_t ei = 0; ei < EN(); ei++) {
elements[ei] = mesh.elements[ei]; elements[ei] = m.elements[ei];
} }
reset(); reset();
}
SimulationMesh::SimulationMesh(VCGTriMesh &triMesh)
{
vcg::tri::Append<VCGEdgeMesh, VCGTriMesh>::MeshCopy(*this, triMesh);
label = triMesh.getLabel();
// eigenEdges = triMesh.getEigenEdges();
// if (eigenEdges.rows() == 0) {
getEdges(eigenEdges);
// }
// eigenVertices = triMesh.getEigenVertices();
// if (eigenVertices.rows() == 0) {
getVertices(eigenVertices);
// }
vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(*this,
std::string(
"Elements"));
nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(*this,
std::string("Nodes"));
initializeNodes();
initializeElements();
reset();
} }
void SimulationMesh::computeElementalProperties() { void SimulationMesh::computeElementalProperties() {
@ -192,15 +218,15 @@ void SimulationMesh::initializeElements() {
2, std::vector<VectorType>(6, VectorType(0, 0, 0))); 2, std::vector<VectorType>(6, VectorType(0, 0, 0)));
element.derivativeT3.resize( element.derivativeT3.resize(
2, std::vector<VectorType>(6, VectorType(0, 0, 0))); 2, std::vector<VectorType>(6, VectorType(0, 0, 0)));
element.derivativeT1_jplus1.resize(6); element.derivativeT1_jplus1.resize(6, VectorType(0, 0, 0));
element.derivativeT1_j.resize(6); element.derivativeT1_j.resize(6, VectorType(0, 0, 0));
element.derivativeT1_jplus1.resize(6); element.derivativeT1_jplus1.resize(6, VectorType(0, 0, 0));
element.derivativeT2_j.resize(6); element.derivativeT2_j.resize(6, VectorType(0, 0, 0));
element.derivativeT2_jplus1.resize(6); element.derivativeT2_jplus1.resize(6, VectorType(0, 0, 0));
element.derivativeT3_j.resize(6); element.derivativeT3_j.resize(6, VectorType(0, 0, 0));
element.derivativeT3_jplus1.resize(6); element.derivativeT3_jplus1.resize(6, VectorType(0, 0, 0));
element.derivativeR_j.resize(6); element.derivativeR_j.resize(6, VectorType(0, 0, 0));
element.derivativeR_jplus1.resize(6); element.derivativeR_jplus1.resize(6, VectorType(0, 0, 0));
} }
updateElementalFrames(); updateElementalFrames();
} }

1
simulationmesh.hpp
View File

@ -28,6 +28,7 @@ public:
SimulationMesh(PatternGeometry &pattern); SimulationMesh(PatternGeometry &pattern);
SimulationMesh(ConstVCGEdgeMesh &edgeMesh); SimulationMesh(ConstVCGEdgeMesh &edgeMesh);
SimulationMesh(SimulationMesh &elementalMesh); SimulationMesh(SimulationMesh &elementalMesh);
SimulationMesh(VCGTriMesh &triMesh);
void updateElementalLengths(); void updateElementalLengths();
void updateIncidentElements(); void updateIncidentElements();

52
topologyenumerator.cpp
View File

@ -5,7 +5,7 @@
#include <numeric> #include <numeric>
#include <unordered_set> #include <unordered_set>
const bool debugIsOn{false}; const bool debugIsOn{true};
const bool savePlyFiles{true}; const bool savePlyFiles{true};
// size_t binomialCoefficient(size_t n, size_t m) { // size_t binomialCoefficient(size_t n, size_t m) {
@ -456,9 +456,8 @@ void TopologyEnumerator::computeValidPatterns(
.string() .string()
+ ".ply"); + ".ply");
} }
} else {
continue; // should be uncommented in order to improve performance
} }
continue; // should be uncommented in order to improve performance
} }
const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold
@ -553,6 +552,50 @@ void TopologyEnumerator::computeValidPatterns(
+ ".ply"); + ".ply");
PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(
patternGeometry); // the marked nodes of hasDanglingEdges are patternGeometry); // the marked nodes of hasDanglingEdges are
std::vector<std::pair<int, PatternGeometry::EdgePointer>> eCC;
vcg::tri::Clean<PatternGeometry>::edgeMeshConnectedComponents(tiledPatternGeometry,
eCC);
vcg::tri::UpdateFlags<PatternGeometry>::EdgeClear(tiledPatternGeometry);
const size_t numberOfCC_edgeBased = eCC.size();
std::sort(eCC.begin(),
eCC.end(),
[](const std::pair<int, PatternGeometry::EdgePointer> &a,
const std::pair<int, PatternGeometry::EdgePointer> &b) {
return a.first > b.first;
});
PatternGeometry::EdgePointer &ep = eCC[0].second;
size_t colorsRegistered = 0;
std::stack<EdgePointer> stack;
stack.push(ep);
while (!stack.empty()) {
EdgePointer ep = stack.top();
stack.pop();
for (int i = 0; i < 2; ++i) {
vcg::edge::VEIterator<PatternGeometry::EdgeType> vei(ep->V(i));
while (!vei.End()) {
if (!vei.E()->IsV()) {
vei.E()->SetV();
stack.push(vei.E());
tiledPatternGeometry
.vert[tiledPatternGeometry.getIndex(vei.V1())]
.C()
= vcg::Color4b::Blue;
tiledPatternGeometry
.vert[tiledPatternGeometry.getIndex(vei.V0())]
.C()
= vcg::Color4b::Blue;
colorsRegistered++;
}
++vei;
}
}
}
assert(colorsRegistered == eCC[0].first);
tiledPatternGeometry.save(std::filesystem::path(moreThanOneCCPath) tiledPatternGeometry.save(std::filesystem::path(moreThanOneCCPath)
.append(patternName + "_tiled") .append(patternName + "_tiled")
.string() .string()
@ -563,8 +606,7 @@ void TopologyEnumerator::computeValidPatterns(
} }
} }
if (hasArticulationPoints /*&& !hasFloatingComponents && if (hasArticulationPoints && !hasFloatingComponents && !tiledPatternHasDanglingEdges) {
!tiledPatternHasDanglingEdges*/) {
statistics.numberOfPatternsWithArticulationPoints++; statistics.numberOfPatternsWithArticulationPoints++;
if (debugIsOn) { if (debugIsOn) {
if (savePlyFiles) { if (savePlyFiles) {

16
trianglepatterngeometry.cpp
View File

@ -112,9 +112,9 @@ PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern) {
vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(tile); vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(tile);
vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(tile); vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(tile);
for (size_t vi = 0; vi < pattern.vn; vi++) { // for (size_t vi = 0; vi < pattern.vn; vi++) {
tile.vert[vi].C() = vcg::Color4b::Blue; // tile.vert[vi].C() = vcg::Color4b::Blue;
} // }
tile.updateEigenEdgeAndVertices(); tile.updateEigenEdgeAndVertices();
return tile; return tile;
@ -816,7 +816,7 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
std::vector<ConstPatternGeometry> &patterns, std::vector<ConstPatternGeometry> &patterns,
const std::vector<int> &connectToNeighborsVi, const std::vector<int> &connectToNeighborsVi,
const VCGPolyMesh &tileInto, const VCGPolyMesh &tileInto,
const std::vector<size_t> &perSurfaceFacePatternIndices, const std::vector<int> &perSurfaceFacePatternIndices,
std::vector<size_t> &tileIntoEdgesToTiledVi, std::vector<size_t> &tileIntoEdgesToTiledVi,
std::vector<std::vector<size_t>> &perPatternIndexToTiledPatternEdgeIndex) std::vector<std::vector<size_t>> &perPatternIndexToTiledPatternEdgeIndex)
{ {
@ -843,6 +843,10 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
assert(vcg::tri::HasFEAdjacency(tileInto)); assert(vcg::tri::HasFEAdjacency(tileInto));
assert(vcg::tri::HasFVAdjacency(tileInto)); assert(vcg::tri::HasFVAdjacency(tileInto));
for (const VCGPolyMesh::FaceType &f : tileInto.face) { for (const VCGPolyMesh::FaceType &f : tileInto.face) {
const int facePatternIndex = perSurfaceFacePatternIndices[tileInto.getIndex(f)];
if (facePatternIndex == -1) {
continue;
}
CoordType centerOfFace(0, 0, 0); CoordType centerOfFace(0, 0, 0);
for (size_t vi = 0; vi < f.VN(); vi++) { for (size_t vi = 0; vi < f.VN(); vi++) {
centerOfFace = centerOfFace + f.cP(vi); centerOfFace = centerOfFace + f.cP(vi);
@ -856,7 +860,6 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
for (int &vi : firstInFanConnectToNeighbor_vi) { for (int &vi : firstInFanConnectToNeighbor_vi) {
vi += pTiledPattern->VN(); vi += pTiledPattern->VN();
} }
const size_t facePatternIndex = perSurfaceFacePatternIndices[tileInto.getIndex(f)];
ConstPatternGeometry &pattern = patterns[facePatternIndex]; ConstPatternGeometry &pattern = patterns[facePatternIndex];
for (size_t vi = 0; vi < f.VN(); vi++) { for (size_t vi = 0; vi < f.VN(); vi++) {
@ -966,6 +969,9 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
tileIntoEdgesToTiledVi.resize(tileInto.EN()); tileIntoEdgesToTiledVi.resize(tileInto.EN());
for (int ei = 0; ei < tileInto.EN(); ei++) { for (int ei = 0; ei < tileInto.EN(); ei++) {
const std::vector<int> interfaceVis = tileIntoEdgeToInterfaceVi[ei]; const std::vector<int> interfaceVis = tileIntoEdgeToInterfaceVi[ei];
if (interfaceVis.empty()) {
continue;
}
assert(interfaceVis.size() == 1 || interfaceVis.size() == 2); assert(interfaceVis.size() == 1 || interfaceVis.size() == 2);
assert( assert(
interfaceVis.size() == 1 interfaceVis.size() == 1

2
trianglepatterngeometry.hpp
View File

@ -120,7 +120,7 @@ private:
std::vector<ConstPatternGeometry> &patterns, std::vector<ConstPatternGeometry> &patterns,
const std::vector<int> &connectToNeighborsVi, const std::vector<int> &connectToNeighborsVi,
const VCGPolyMesh &tileInto, const VCGPolyMesh &tileInto,
const std::vector<size_t> &perSurfaceFacePatternIndices, const std::vector<int> &perSurfaceFacePatternIndices,
std::vector<size_t> &tileIntoEdgesToTiledVi, std::vector<size_t> &tileIntoEdgesToTiledVi,
std::vector<std::vector<size_t>> &perPatternIndexTiledPatternEdgeIndex); std::vector<std::vector<size_t>> &perPatternIndexTiledPatternEdgeIndex);
}; };

2
utilities.hpp
View File

@ -315,7 +315,7 @@ size_t computeHashUnordered(const std::vector<T> &v)
return hash; return hash;
} }
inline size_t computeHashOrdered(const std::vector<size_t> &v) inline size_t computeHashOrdered(const std::vector<int> &v)
{ {
std::string elementsString; std::string elementsString;
for (const auto &el : v) { for (const auto &el : v) {

35
vcgtrimesh.cpp
View File

@ -8,27 +8,31 @@
bool VCGTriMesh::load(const std::string &filename) { bool VCGTriMesh::load(const std::string &filename) {
assert(std::filesystem::exists(filename)); assert(std::filesystem::exists(filename));
// unsigned int mask = 0; unsigned int mask = 0;
// mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTCOORD; mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTCOORD;
// mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTNORMAL; mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTNORMAL;
// mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTCOLOR; mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTCOLOR;
// mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_EDGEINDEX; mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_EDGEINDEX;
// mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_FACEINDEX; mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_FACEINDEX;
// if (nanoply::NanoPlyWrapper<VCGTriMesh>::LoadModel( if (nanoply::NanoPlyWrapper<VCGTriMesh>::LoadModel(
// std::filesystem::absolute(filename).string().c_str(), *this, mask) != 0) { std::filesystem::absolute(filename).string().c_str(), *this, mask)
if (vcg::tri::io::Importer<VCGTriMesh>::Open(*this,
std::filesystem::absolute(filename).string().c_str())
!= 0) { != 0) {
// if (vcg::tri::io::Importer<VCGTriMesh>::Open(*this,
// std::filesystem::absolute(filename).string().c_str())
// != 0) {
std::cout << "Could not load tri mesh" << std::endl; std::cout << "Could not load tri mesh" << std::endl;
return false; return false;
} }
vcg::tri::UpdateTopology<VCGTriMesh>::FaceFace(*this);
vcg::tri::UpdateTopology<VCGTriMesh>::VertexFace(*this);
vcg::tri::UpdateNormal<VCGTriMesh>::PerVertexNormalized(*this);
label = std::filesystem::path(filename).stem().string(); vcg::tri::UpdateTopology<VCGTriMesh>::AllocateEdge(*this);
return true; vcg::tri::UpdateTopology<VCGTriMesh>::FaceFace(*this);
vcg::tri::UpdateTopology<VCGTriMesh>::VertexFace(*this);
vcg::tri::UpdateTopology<VCGTriMesh>::VertexEdge(*this);
vcg::tri::UpdateNormal<VCGTriMesh>::PerVertexNormalized(*this);
label = std::filesystem::path(filename).stem().string();
return true;
} }
Eigen::MatrixX3d VCGTriMesh::getVertices() const Eigen::MatrixX3d VCGTriMesh::getVertices() const
@ -87,6 +91,7 @@ bool VCGTriMesh::save(const std::string plyFilename)
} }
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
polyscope::SurfaceMesh *VCGTriMesh::registerForDrawing( polyscope::SurfaceMesh *VCGTriMesh::registerForDrawing(
const std::optional<std::array<double, 3>> &desiredColor, const bool &shouldEnable) const const std::optional<std::array<double, 3>> &desiredColor, const bool &shouldEnable) const
{ {

8
vcgtrimesh.hpp
View File

@ -21,14 +21,15 @@ class VCGTriMeshVertex : public vcg::Vertex<VCGTriMeshUsedTypes,
vcg::vertex::BitFlags, vcg::vertex::BitFlags,
vcg::vertex::Color4b, vcg::vertex::Color4b,
vcg::vertex::VFAdj, vcg::vertex::VFAdj,
vcg::vertex::Qualityd> vcg::vertex::Qualityd,
vcg::vertex::VEAdj>
{}; {};
class VCGTriMeshFace class VCGTriMeshFace
: public vcg::Face<VCGTriMeshUsedTypes, vcg::face::FFAdj, vcg::face::VFAdj, : public vcg::Face<VCGTriMeshUsedTypes, vcg::face::FFAdj, vcg::face::VFAdj,
vcg::face::VertexRef, vcg::face::BitFlags, vcg::face::VertexRef, vcg::face::BitFlags,
vcg::face::Normal3d> {}; vcg::face::Normal3d> {};
class VCGTriMeshEdge class VCGTriMeshEdge : public vcg::Edge<VCGTriMeshUsedTypes, vcg::edge::VertexRef, vcg::edge::VEAdj>
: public vcg::Edge<VCGTriMeshUsedTypes, vcg::edge::VertexRef> {}; {};
class VCGTriMesh : public vcg::tri::TriMesh<std::vector<VCGTriMeshVertex>, class VCGTriMesh : public vcg::tri::TriMesh<std::vector<VCGTriMeshVertex>,
std::vector<VCGTriMeshFace>, std::vector<VCGTriMeshFace>,
@ -52,6 +53,7 @@ public:
const bool &shouldEnable = true) const; const bool &shouldEnable = true) const;
#endif #endif
Eigen::MatrixX2i getEdges() const; Eigen::MatrixX2i getEdges() const;
void updateEigenEdgeAndVertices();
}; };
#endif // VCGTRIMESH_HPP #endif // VCGTRIMESH_HPP