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

Refacoring. Drm has res forces norm threshold 1e-3

This commit is contained in:
iasonmanolas 2021-06-24 10:02:20 +03:00
parent 33774c546c
commit 3e3d32399b
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.xi = 0.9969;
settings.maxDRMIterations = 0;
formFinder.totalResidualForcesNormThreshold = 1;
formFinder.mSettings.totalResidualForcesNormThreshold = 1;
settings.shouldDraw = false;
settings.beVerbose = false;
// settings.shouldCreatePlots = true;
settings.beVerbose = true;
settings.shouldCreatePlots = true;
SimulationResults simpleBeam_simulationResults
= formFinder.executeSimulation(std::make_shared<SimulationJob>(beamSimulationJob), settings);
simpleBeam_simulationResults.save();
@ -227,11 +227,10 @@ void DRMSimulationModel::reset()
plotYValues.clear();
plotHandle.reset();
checkedForMaximumMoment = false;
mSettings.shouldUseTranslationalKineticEnergyThreshold = false;
externalMomentsNorm = 0;
mSettings.drawingStep = 1;
Dt = mSettings.Dtini;
numOfDampings = 0;
shouldTemporarilyDampForces = false;
}
VectorType DRMSimulationModel::computeDisplacementDifferenceDerivative(
@ -539,6 +538,7 @@ double DRMSimulationModel::computeTheta3(const EdgeType &e, const VertexType &v)
// Should be refactored in the future
double theta3;
const bool shouldBreak = mCurrentSimulationStep == 12970;
if (&e == node.referenceElement) {
// Use nR as theta3 only for the first star edge
return nR;
@ -794,7 +794,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
// omp_lock_t writelock;
// omp_init_lock(&writelock);
#ifdef ENABLE_OPENMP
#pragma omp parallel for schedule(static) num_threads(5)
#pragma omp parallel for schedule(static) num_threads(8)
#endif
for (int ei = 0; ei < pMesh->EN(); 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 theta2_j = ef.t2.dot(t3CrossN_j);
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_jplus1 = computeTheta3(e, ev_jplus1);
// element.rotationalDisplacements_j.theta1 = theta1_j;
@ -837,7 +838,6 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
const bool isDofConstrainedFor_ev = fixedVertices.contains(edgeNode.vi)
&& fixedVertices.at(edgeNode.vi).contains(dofi);
const bool shouldBreak = edgeNode.vi == 0 && dofi == 5;
if (!isDofConstrainedFor_ev) {
DifferentiateWithRespectTo dui{ev, dofi};
// Axial force computation
@ -976,8 +976,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
for (size_t vi = 0; vi < pMesh->VN(); vi++) {
Node::Forces &force = pMesh->nodes[vi].force;
const Vector6d &nodeResidualForce = force.residual;
// sumOfResidualForces = sumOfResidualForces + nodeResidualForce;
const double residualForceNorm = nodeResidualForce.norm();
sumOfResidualForces = sumOfResidualForces + nodeResidualForce;
// const double residualForceNorm = nodeResidualForce.norm();
const double residualForceNorm = nodeResidualForce.getTranslation().norm();
const bool shouldTrigger = std::isnan(residualForceNorm);
totalResidualForcesNorm += residualForceNorm;
}
@ -1145,7 +1146,11 @@ DRMSimulationModel::DRMSimulationModel() {}
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) {
const size_t vi = pMesh->getIndex(v);
double translationalSumSk = 0;
@ -1208,6 +1213,14 @@ void DRMSimulationModel::updateNodalAccelerations()
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);
Node &node = pMesh->nodes[v];
node.velocity = node.velocity + node.acceleration * Dt;
// if (vi == 10) {
// std::cout << "Velocity:" << node.velocity[0] << " " << node.velocity[1] << " "
// << node.velocity[2] << std::endl;
// }
}
updateKineticEnergy();
}
void DRMSimulationModel::updateNodalDisplacements()
{
const bool shouldCapDisplacements = mSettings.displacementCap.has_value();
for (VertexType &v : pMesh->vert) {
Node &node = pMesh->nodes[v];
node.displacements = node.displacements + node.velocity * Dt;
// if (mSettings.beVerbose) {
// 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;
// }
Vector6d stepDisplacement = node.velocity * Dt;
if (shouldCapDisplacements
&& stepDisplacement.getTranslation().norm() > mSettings.displacementCap) {
stepDisplacement = stepDisplacement
* (*mSettings.displacementCap
/ stepDisplacement.getTranslation().norm());
std::cout << "Displacement capped" << std::endl;
}
node.displacements = node.displacements + stepDisplacement;
// 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(
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 (!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) {
// std::cout << "AFTER:" << mesh->vert[1].N()[0] << " " <<
// mesh->vert[1].N()[1]
@ -1345,6 +1372,7 @@ void DRMSimulationModel::applyDisplacements(
for (VertexType &v : pMesh->vert) {
updateNodePosition(v, fixedVertices);
updateNodeNormal(v, fixedVertices);
updateNodeNr(v);
}
updateElementalFrames();
if (mSettings.shouldDraw) {
@ -1427,8 +1455,8 @@ void DRMSimulationModel::updateKineticEnergy()
const double nodeRotationalKineticEnergy
= 0.5
* (node.mass.rotationalJ * pow(node.velocity[3], 2)
+ +node.mass.rotationalI3 * pow(node.velocity[4], 2)
+ +node.mass.rotationalI2 * pow(node.velocity[5], 2));
+ node.mass.rotationalI3 * pow(node.velocity[4], 2)
+ node.mass.rotationalI2 * pow(node.velocity[5], 2));
node.kineticEnergy += nodeTranslationalKineticEnergy /*+ nodeRotationalKineticEnergy*/;
assert(node.kineticEnergy < 1e15);
@ -1442,12 +1470,13 @@ void DRMSimulationModel::updateKineticEnergy()
void DRMSimulationModel::resetVelocities()
{
for (const VertexType &v : pMesh->vert) {
pMesh->nodes[v].velocity = 0;
// mesh->nodes[v].force.residual * 0.5 * Dt /
// mesh->nodes[v].mass; // NOTE: Do I reset the previous
// velocity?
// reset
// current to 0 or this?
pMesh->nodes[v].velocity =
// pMesh->nodes[v].acceleration
// * Dt; // NOTE: Do I reset the previous
// velocity?
// reset
// current to 0 or this?
0;
}
updateKineticEnergy();
}
@ -1584,6 +1613,9 @@ void DRMSimulationModel::printCurrentState() const
* 1e-6 / (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;
}
void DRMSimulationModel::printDebugInfo() const
@ -1660,12 +1692,13 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
accelerationX[pMesh->getIndex(v)] = pMesh->nodes[v].acceleration[0];
}
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("External force", externalForces)->setEnabled(false);
meshPolyscopeHandle->addNodeVectorQuantity("External Moments", externalMoments)->setEnabled(true);
meshPolyscopeHandle->addNodeScalarQuantity("Internal force norm", internalForcesNorm)
->setEnabled(true);
meshPolyscopeHandle->setRadius(0.002);
// polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeVectorQuantity("Internal Axial force", internalAxialForces)
// ->setEnabled(false);
@ -1735,8 +1768,9 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
// instead of full width. Must have
// matching PopItemWidth() below.
ImGui::InputInt("Simulation drawing step",
&mSettings.drawingStep); // set a int variable
ImGui::InputInt("Simulation debug step",
&mSettings.debugModeStep); // set a int variable
mSettings.drawingStep = mSettings.debugModeStep;
ImGui::Checkbox("Enable drawing",
&mSettings.shouldDraw); // set a int variable
ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep);
@ -1766,30 +1800,183 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
}
#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,
const Settings &settings)
const Settings &settings,
const SimulationResults &solutionGuess)
{
assert(pJob->pMesh.operator bool());
auto t1 = std::chrono::high_resolution_clock::now();
reset();
mSettings = settings;
double totalExternalForcesNorm = 0;
for (const std::pair<VertexIndex, Vector6d> &externalForce : pJob->nodalExternalForces) {
totalExternalForcesNorm += externalForce.second.norm();
}
reset();
// double totalExternalForcesNorm = 0;
// for (const std::pair<VertexIndex, Vector6d> &externalForce : pJob->nodalExternalForces) {
// totalExternalForcesNorm += externalForce.second.norm();
// }
if (!pJob->nodalExternalForces.empty()) {
totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination
* totalExternalForcesNorm;
} else {
totalResidualForcesNormThreshold = 1;
std::cout << "No forces setted default residual forces norm threshold" << std::endl;
}
// if (!pJob->nodalExternalForces.empty()) {
// // totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination
// // * totalExternalForcesNorm;
// } else {
// totalResidualForcesNormThreshold = 1;
// std::cout << "No forces setted default residual forces norm threshold" << std::endl;
// }
if (mSettings.beVerbose) {
std::cout << "totalResidualForcesNormThreshold:" << totalResidualForcesNormThreshold
<< std::endl;
std::cout << "totalResidualForcesNormThreshold:"
<< mSettings.totalResidualForcesNormThreshold << std::endl;
}
history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
// if (!pJob->nodalExternalForces.empty()) {
// double externalForcesNorm = 0;
@ -1835,141 +2022,170 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
}
#endif
updateElementalFrames();
updateNodalMasses();
if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) {
shouldApplyInitialDistortion = true;
}
updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
const size_t movingAverageSampleSize = 50;
std::queue<double> residualForcesMovingAverageHistorySample;
if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) {
shouldApplyInitialDistortion = true;
}
updateElementalFrames();
if (!solutionGuess.displacements.empty()) {
//#ifdef POLYSCOPE_DEFINED
// if (mSettings.shouldDraw || mSettings.isDebugMode) {
// solutionGuess.registerForDrawing();
// polyscope::show();
// solutionGuess.unregister();
// }
//#endif
double residualForcesMovingAverageDerivativeMax = 0;
while (settings.maxDRMIterations == 0 || mCurrentSimulationStep < settings.maxDRMIterations) {
// while (true) {
updateNormalDerivatives();
updateT1Derivatives();
updateRDerivatives();
updateT2Derivatives();
updateT3Derivatives();
updateResidualForcesOnTheFly(constrainedVertices);
applySolutionGuess(solutionGuess, pJob);
shouldTemporarilyDampForces = true;
}
// TODO: write parallel function for updating positions
// TODO: Make a single function out of updateResidualForcesOnTheFly
// updatePositionsOnTheFly
// updatePositionsOnTheFly(constrainedVertices);
updateNodalMasses();
updateNodalAccelerations();
updateNodalVelocities();
updateNodalDisplacements();
applyDisplacements(constrainedVertices);
if (!pJob->nodalForcedDisplacements.empty()) {
applyForcedDisplacements(
updateNodalMasses();
updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
const size_t movingAverageSampleSize = 200;
std::queue<double> residualForcesMovingAverageHistorySample;
std::vector<double> percentageOfConvergence;
pJob->nodalForcedDisplacements);
}
// if (!pJob->nodalForcedNormals.empty()) {
// applyForcedNormals(pJob->nodalForcedNormals);
// }
updateElementalLengths();
mCurrentSimulationStep++;
double residualForcesMovingAverageDerivativeMax = 0;
while (settings.maxDRMIterations == 0 || mCurrentSimulationStep < settings.maxDRMIterations) {
// while (true) {
updateNormalDerivatives();
updateT1Derivatives();
updateRDerivatives();
updateT2Derivatives();
updateT3Derivatives();
const bool shouldBreak = mCurrentSimulationStep == 12970;
updateResidualForcesOnTheFly(constrainedVertices);
double sumOfDisplacementsNorm = 0;
for (size_t vi = 0; vi < pMesh->VN(); vi++) {
sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm();
}
sumOfDisplacementsNorm /= pMesh->bbox.Diag();
pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm;
// pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm
// + mCurrentSimulationStep
// * pMesh->residualForcesMovingAverage)
// / (1 + mCurrentSimulationStep);
pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm
/ movingAverageSampleSize;
residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage);
if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) {
const double firstElementValue = residualForcesMovingAverageHistorySample.front();
pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize;
residualForcesMovingAverageHistorySample.pop();
// TODO: write parallel function for updating positions
// TODO: Make a single function out of updateResidualForcesOnTheFly
// updatePositionsOnTheFly
// updatePositionsOnTheFly(constrainedVertices);
updateNodalMasses();
updateNodalAccelerations();
updateNodalVelocities();
updateNodalDisplacements();
applyDisplacements(constrainedVertices);
if (!pJob->nodalForcedDisplacements.empty()) {
applyForcedDisplacements(
pMesh->residualForcesMovingAverageDerivativeNorm
= std::abs((residualForcesMovingAverageHistorySample.back()
- residualForcesMovingAverageHistorySample.front()))
/ (movingAverageSampleSize);
residualForcesMovingAverageDerivativeMax
= std::max(pMesh->residualForcesMovingAverageDerivativeNorm,
residualForcesMovingAverageDerivativeMax);
pMesh->residualForcesMovingAverageDerivativeNorm
/= residualForcesMovingAverageDerivativeMax;
// std::cout << "Normalized derivative:"
// << residualForcesMovingAverageDerivativeNorm
// << std::endl;
}
pJob->nodalForcedDisplacements);
}
// if (!pJob->nodalForcedNormals.empty()) {
// applyForcedNormals(pJob->nodalForcedNormals);
// }
updateElementalLengths();
mCurrentSimulationStep++;
// pMesh->previousTotalPotentialEnergykN =
// pMesh->currentTotalPotentialEnergykN;
// pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000;
// timePerNodePerIterationHistor.push_back(timePerNodePerIteration);
if (mSettings.beVerbose && mSettings.isDebugMode
&& mCurrentSimulationStep % mSettings.debugModeStep == 0) {
printCurrentState();
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver",
// movingAverages);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver deriv",
// movingAveragesDerivatives);
// draw();
std::cout << "Mov aver deriv:" << pMesh->residualForcesMovingAverageDerivativeNorm
<< std::endl;
// SimulationResulnodalForcedDisplacementstsReporter::createPlot("Number of Steps",
// "Time/(#nodes*#iterations)",
// timePerNodePerIterationHistory);
}
double sumOfDisplacementsNorm = 0;
for (size_t vi = 0; vi < pMesh->VN(); vi++) {
sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm();
}
sumOfDisplacementsNorm /= pMesh->bbox.Diag();
pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm;
// pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm
// + mCurrentSimulationStep
// * pMesh->residualForcesMovingAverage)
// / (1 + mCurrentSimulationStep);
pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm
/ movingAverageSampleSize;
residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage);
if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) {
const double firstElementValue = residualForcesMovingAverageHistorySample.front();
pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize;
residualForcesMovingAverageHistorySample.pop();
if ((mSettings.shouldCreatePlots || mSettings.isDebugMode) && mCurrentSimulationStep != 0) {
history.stepPulse(*pMesh);
}
pMesh->residualForcesMovingAverageDerivativeNorm
= std::abs((residualForcesMovingAverageHistorySample.back()
- residualForcesMovingAverageHistorySample.front()))
/ (movingAverageSampleSize);
residualForcesMovingAverageDerivativeMax
= std::max(pMesh->residualForcesMovingAverageDerivativeNorm,
residualForcesMovingAverageDerivativeMax);
pMesh->residualForcesMovingAverageDerivativeNorm
/= residualForcesMovingAverageDerivativeMax;
// std::cout << "Normalized derivative:"
// << residualForcesMovingAverageDerivativeNorm
// << std::endl;
}
if (std::isnan(pMesh->currentTotalKineticEnergy)) {
if (mSettings.beVerbose) {
std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl;
}
break;
}
// pMesh->previousTotalPotentialEnergykN =
// pMesh->currentTotalPotentialEnergykN;
// pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000;
// timePerNodePerIterationHistor.push_back(timePerNodePerIteration);
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
if ((mSettings.shouldDraw && mCurrentSimulationStep % mSettings.drawingStep == 0)
|| (mSettings.isDebugMode && mCurrentSimulationStep % mSettings.debugModeStep == 0)) /* &&
if (mSettings.shouldDraw && mSettings.isDebugMode
&& mCurrentSimulationStep % mSettings.drawingStep == 0) /* &&
currentSimulationStep > maxDRMIterations*/
{
// std::string saveTo = std::filesystem::current_path()
// .append("Debugging_files")
// .append("Screenshots")
// .string();
draw();
// yValues = history.kineticEnergy;
// plotHandle = matplot::scatter(xPlot, yValues);
// label = "Log of Kinetic energy";
// plotHandle->legend_string(label);
{
// std::string saveTo = std::filesystem::current_path()
// .append("Debugging_files")
// .append("Screenshots")
// .string();
draw();
// yValues = history.kineticEnergy;
// plotHandle = matplot::scatter(xPlot, yValues);
// label = "Log of Kinetic energy";
// plotHandle->legend_string(label);
// shouldUseKineticEnergyThreshold = true;
}
// shouldUseKineticEnergyThreshold = true;
}
#endif
if (mSettings.isDebugMode && mCurrentSimulationStep % mSettings.debugModeStep == 0) {
if (mSettings.shouldCreatePlots && mSettings.isDebugMode
&& mCurrentSimulationStep % mSettings.debugModeStep == 0) {
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver deriv",
// 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",
// "Log of Kinetic energy",
// history.kineticEnergy);
// "Residual Forces mov aver",
// history.residualForcesMovingAverage);
SimulationResultsReporter::createPlot("Number of Steps",
"Log of Residual Forces",
history.logResidualForces,
{},
history.redMarks);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Kinetic energy",
// history.kineticEnergy);
// SimulationResultsReporter reporter;
// reporter.reportHistory(history, "Results");
// SimulationResultsReporter::createPlot("Number of Steps",
// "Percentage of convergence",
// percentageOfConvergence,
// {},
// history.redMarks);
}
// t = t + Dt;
// std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy
@ -1990,16 +2206,19 @@ mesh->currentTotalPotentialEnergykN*/
= mSettings.shouldUseTranslationalKineticEnergyThreshold
&& pMesh->currentTotalTranslationalKineticEnergy
< mSettings.totalTranslationalKineticEnergyThreshold
&& mCurrentSimulationStep > 20 && numOfDampings > 2;
&& mCurrentSimulationStep > 20 && numOfDampings > 0;
const bool fullfillsResidualForcesNormTerminationCriterion
= pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold;
= pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold;
const bool fullfillsMovingAverageNormTerminationCriterion
= pMesh->residualForcesMovingAverage
< mSettings.residualForcesMovingAverageNormThreshold;
/*pMesh->residualForcesMovingAverage < totalResidualForcesNormThreshold*/
const bool fullfillsMovingAverageDerivativeNormTerminationCriterion
= pMesh->residualForcesMovingAverageDerivativeNorm < 1e-4;
const bool shouldTerminate
= fullfillsKineticEnergyTerminationCriterion
// || fullfillsMovingAverageNormTerminationCriterion
// || fullfillsMovingAverageDerivativeNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion;
if (shouldTerminate) {
if (mSettings.beVerbose) {
@ -2021,34 +2240,59 @@ mesh->currentTotalPotentialEnergykN*/
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);
// }
// updateElementalLengths();
// if (mSettings.beVerbose) {
// printCurrentState();
// }
// 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();
Dt = Dt * mSettings.xi;
++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);
if (mCurrentSimulationStep == settings.maxDRMIterations &&
mCurrentSimulationStep != 0) {
std::cout << "Did not reach equilibrium before reaching the maximum number "
"of DRM steps ("
<< settings.maxDRMIterations << "). Breaking simulation"
<< std::endl;
results.converged = false;
if (mCurrentSimulationStep == settings.maxDRMIterations && mCurrentSimulationStep != 0) {
std::cout << "Did not reach equilibrium before reaching the maximum number "
"of DRM steps ("
<< settings.maxDRMIterations << "). Breaking simulation" << std::endl;
results.converged = false;
} else if (std::isnan(pMesh->currentTotalKineticEnergy)) {
std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl;
results.converged = false;
@ -2064,7 +2308,6 @@ mesh->currentTotalPotentialEnergykN*/
<< "s" << std::endl;
std::cout << "Number of dampings:" << numOfDampings << endl;
std::cout << "Number of steps:" << mCurrentSimulationStep << endl;
results.converged = true;
}
// mesh.printVertexCoordinates(mesh.VN() / 2);
#ifdef POLYSCOPE_DEFINED
@ -2073,46 +2316,93 @@ mesh->currentTotalPotentialEnergykN*/
}
#endif
if (results.converged) {
results.debug_drmDisplacements = results.displacements;
results.internalPotentialEnergy = computeTotalInternalPotentialEnergy();
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++) {
const Node &node = pMesh->nodes[vi];
const Eigen::Vector3d initialNormal = node.initialNormal.ToEigenVector<Eigen::Vector3d>();
const Eigen::Vector3d deformedNormal
const Eigen::Vector3d nInitial_eigen = node.initialNormal.ToEigenVector<Eigen::Vector3d>();
const Eigen::Vector3d nDeformed_eigen
= pMesh->vert[vi].cN().ToEigenVector<Eigen::Vector3d>();
Eigen::Quaternion<double> q;
q.setFromTwoVectors(initialNormal, deformedNormal);
Eigen::Quaternion<double> q_nr;
q_nr = Eigen::AngleAxis<double>(pMesh->nodes[vi].nR, deformedNormal);
Eigen::Quaternion<double> q_normal;
q_normal.setFromTwoVectors(nInitial_eigen, nDeformed_eigen);
Eigen::Quaternion<double> q_nr_nDeformed;
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;
const VectorType &n_deformed = pMesh->vert[vi].cN();
const VectorType &nDeformed = pMesh->vert[vi].cN();
const VectorType referenceF1_deformed = (referenceT1_deformed
- (n_deformed
* (referenceT1_deformed * n_deformed)))
- (node.initialNormal
* (referenceT1_deformed * node.initialNormal)))
.Normalize();
const VectorType referenceT1_initial
= computeT1Vector(pMesh->nodes[node.referenceElement->cV(0)].initialLocation,
pMesh->nodes[node.referenceElement->cV(1)].initialLocation);
const VectorType &n_initial = node.initialNormal;
const VectorType referenceF1_initial
= (referenceT1_initial - (n_deformed * (referenceT1_initial * n_deformed))).Normalize();
Eigen::Quaternion<double> q_f1; //nr is with respect to f1
q_f1.setFromTwoVectors(referenceF1_initial.ToEigenVector<Eigen::Vector3d>(),
referenceF1_deformed.ToEigenVector<Eigen::Vector3d>());
Eigen::Quaternion<double> q_t1;
q_t1.setFromTwoVectors(referenceT1_initial.ToEigenVector<Eigen::Vector3d>(),
referenceT1_deformed.ToEigenVector<Eigen::Vector3d>());
// const VectorType &n_initial = node.initialNormal;
const VectorType referenceF1_initial = (referenceT1_initial
- (node.initialNormal
* (referenceT1_initial * node.initialNormal)))
.Normalize();
Eigen::Quaternion<double> q_f1_nInit; //nr is with respect to f1
q_f1_nInit.setFromTwoVectors(referenceF1_initial.ToEigenVector<Eigen::Vector3d>(),
referenceF1_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
const Eigen::Vector3d eulerAngles
= results.rotationalDisplacementQuaternion[vi].toRotationMatrix().eulerAngles(0, 1, 2);
results.displacements[vi][3] = eulerAngles[0];
results.displacements[vi][4] = eulerAngles[1];
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
{
bool isDebugMode{false};
int debugModeStep{10000};
int debugModeStep{100000};
bool shouldDraw{false};
bool beVerbose{false};
bool shouldCreatePlots{false};
@ -40,79 +40,78 @@ public:
int maxDRMIterations{0};
bool shouldUseTranslationalKineticEnergyThreshold{false};
int gradualForcedDisplacementSteps{100};
double gamma{0.8};
std::optional<double> displacementCap;
double totalResidualForcesNormThreshold{1e-3};
Settings() {}
};
double totalResidualForcesNormThreshold{1};
private:
Settings mSettings;
double Dt{mSettings.Dtini};
bool checkedForMaximumMoment{false};
double externalMomentsNorm{0};
size_t mCurrentSimulationStep{0};
matplot::line_handle plotHandle;
std::vector<double> plotYValues;
size_t numOfDampings{0};
private:
Settings mSettings;
double Dt{mSettings.Dtini};
bool checkedForMaximumMoment{false};
bool shouldTemporarilyDampForces{false};
bool shouldTemporarilyAmplifyForces{true};
double externalMomentsNorm{0};
size_t mCurrentSimulationStep{0};
matplot::line_handle plotHandle;
std::vector<double> plotYValues;
size_t numOfDampings{0};
const std::string meshPolyscopeLabel{"Simulation mesh"};
std::unique_ptr<SimulationMesh> pMesh;
std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
constrainedVertices;
SimulationHistory history;
// Eigen::Tensor<double, 4> theta3Derivatives;
// std::unordered_map<MyKeyType, double, key_hash> theta3Derivatives;
bool shouldApplyInitialDistortion = false;
std::unordered_set<VertexIndex> rigidSupports;
const std::string meshPolyscopeLabel{"Simulation mesh"};
std::unique_ptr<SimulationMesh> pMesh;
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> constrainedVertices;
SimulationHistory history;
// Eigen::Tensor<double, 4> theta3Derivatives;
// std::unordered_map<MyKeyType, double, key_hash> theta3Derivatives;
bool shouldApplyInitialDistortion = false;
std::unordered_set<VertexIndex> rigidSupports;
void reset();
void updateNodalInternalForces(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices);
void updateNodalExternalForces(
const std::unordered_map<VertexIndex, Vector6d> &nodalForces,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices);
void updateResidualForces();
void updateRotationalDisplacements();
void updateElementalLengths();
void reset();
void updateNodalInternalForces(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void updateNodalExternalForces(
const std::unordered_map<VertexIndex, Vector6d> &nodalForces,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void updateResidualForces();
void updateRotationalDisplacements();
void updateElementalLengths();
void updateNodalMasses();
void updateNodalMasses();
void updateNodalAccelerations();
void updateNodalAccelerations();
void updateNodalVelocities();
void updateNodalVelocities();
void updateNodalDisplacements();
void updateNodalDisplacements();
void applyForcedDisplacements(
const std::unordered_map<VertexIndex, Eigen::Vector3d> &nodalForcedDisplacements);
void applyForcedDisplacements(
const std::unordered_map<VertexIndex, Eigen::Vector3d> &nodalForcedDisplacements);
Vector6d computeElementTorsionalForce(
const Element &element, const Vector6d &displacementDifference,
const std::unordered_set<DoFType> &constrainedDof);
Vector6d computeElementTorsionalForce(const Element &element,
const Vector6d &displacementDifference,
const std::unordered_set<DoFType> &constrainedDof);
// BeamFormFinder::Vector6d computeElementFirstBendingForce(
// const Element &element, const Vector6d &displacementDifference,
// const std::unordered_set<gsl::index> &constrainedDof);
// BeamFormFinder::Vector6d computeElementFirstBendingForce(
// const Element &element, const Vector6d &displacementDifference,
// const std::unordered_set<gsl::index> &constrainedDof);
// BeamFormFinder::Vector6d computeElementSecondBendingForce(
// const Element &element, const Vector6d &displacementDifference,
// const std::unordered_set<gsl::index> &constrainedDof);
// BeamFormFinder::Vector6d computeElementSecondBendingForce(
// const Element &element, const Vector6d &displacementDifference,
// 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(
VertexType &v,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices);
void updateNodePosition(
VertexType &v,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void applyDisplacements(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices);
void applyDisplacements(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
#ifdef POLYSCOPE_DEFINED
void draw(const string &screenshotsFolder= {});
@ -206,11 +205,16 @@ public:
double computeTotalInternalPotentialEnergy();
void applySolutionGuess(const SimulationResults &solutionGuess,
const std::shared_ptr<SimulationJob> &pJob);
void updateNodeNr(VertexType &v);
public:
DRMSimulationModel();
SimulationResults
executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings = Settings());
SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings = Settings(),
const SimulationResults &solutionGuess = SimulationResults());
static void runUnitTests();
};

1
edgemesh.cpp
View File

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

7
polymesh.hpp
View File

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

64
reducedmodeloptimizer_structs.hpp
View File

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

48
simulation_structs.hpp
View File

@ -35,7 +35,7 @@ struct SimulationHistory {
size_t numberOfSteps{0};
std::string label;
std::vector<double> residualForcesLog;
std::vector<double> logResidualForces;
std::vector<double> kineticEnergy;
std::vector<double> potentialEnergies;
std::vector<size_t> redMarks;
@ -48,19 +48,20 @@ struct SimulationHistory {
void markGreen(const size_t &stepNumber) { greenMarks.push_back(stepNumber); }
void stepPulse(const SimulationMesh &mesh) {
kineticEnergy.push_back(log(mesh.currentTotalKineticEnergy));
// potentialEnergy.push_back(mesh.totalPotentialEnergykN);
residualForcesLog.push_back(std::log(mesh.totalResidualForcesNorm));
residualForcesMovingAverage.push_back(mesh.residualForcesMovingAverage);
sumOfNormalizedDisplacementNorms.push_back(mesh.sumOfNormalizedDisplacementNorms);
// residualForcesMovingAverageDerivativesLog.push_back(
// std::log(mesh.residualForcesMovingAverageDerivativeNorm));
void stepPulse(const SimulationMesh &mesh)
{
kineticEnergy.push_back(std::log10(mesh.currentTotalKineticEnergy));
// potentialEnergy.push_back(mesh.totalPotentialEnergykN);
logResidualForces.push_back(std::log10(mesh.totalResidualForcesNorm));
residualForcesMovingAverage.push_back(mesh.residualForcesMovingAverage);
sumOfNormalizedDisplacementNorms.push_back(mesh.sumOfNormalizedDisplacementNorms);
// residualForcesMovingAverageDerivativesLog.push_back(
// std::log(mesh.residualForcesMovingAverageDerivativeNorm));
}
void clear()
{
residualForcesLog.clear();
logResidualForces.clear();
kineticEnergy.clear();
potentialEnergies.clear();
residualForcesMovingAverage.clear();
@ -119,6 +120,11 @@ public:
SimulationJob() {}
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,
SimulationJob &simulationJobMapped) const
{
@ -421,14 +427,18 @@ public:
->setEnabled(shouldEnable);
}
// per node external moments
bool hasExternalMoments = false;
std::vector<std::array<double, 3>> externalMoments(pMesh->VN());
for (const auto &forcePair : nodalExternalForces) {
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]};
}
if (!externalMoments.empty()) {
if (hasExternalMoments) {
polyscope::getCurveNetwork(meshLabel)
->addNodeVectorQuantity("External moment_" + label, externalMoments)
->setEnabled(shouldEnable);
@ -439,8 +449,12 @@ public:
if (polyscope::getCurveNetwork(meshLabel) == nullptr) {
return;
}
polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force_" + label);
polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label);
if (!nodalExternalForces.empty()) {
polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force_" + label);
}
if (!constrainedVertices.empty()) {
polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label);
}
}
#endif // POLYSCOPE_DEFINED
};
@ -449,9 +463,13 @@ struct SimulationResults
/*TODO: remove rotationalDisplacementQuaternion since the last three components of the displacments
* vector contains the same info using euler angles
*/
bool converged{true};
bool converged{false};
std::shared_ptr<SimulationJob> job;
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<Eigen::Quaternion<double>> rotationalDisplacementQuaternion; //per vertex
double internalPotentialEnergy{0};

85
simulationhistoryplotter.hpp
View File

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

68
simulationmesh.cpp
View File

@ -40,20 +40,46 @@ SimulationMesh::SimulationMesh(PatternGeometry &pattern) {
initializeElements();
}
SimulationMesh::SimulationMesh(SimulationMesh &mesh) {
vcg::tri::MeshAssert<SimulationMesh>::VertexNormalNormalized(mesh);
VCGEdgeMesh::copy(mesh);
SimulationMesh::SimulationMesh(SimulationMesh &m)
{
vcg::tri::MeshAssert<SimulationMesh>::VertexNormalNormalized(m);
VCGEdgeMesh::copy(m);
elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(
*this, std::string("Elements"));
nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(
*this, std::string("Nodes"));
initializeNodes();
elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(*this,
std::string(
"Elements"));
nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(*this,
std::string("Nodes"));
initializeNodes();
for (size_t ei = 0; ei < EN(); ei++) {
elements[ei] = mesh.elements[ei];
}
reset();
for (size_t ei = 0; ei < EN(); ei++) {
elements[ei] = m.elements[ei];
}
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() {
@ -192,15 +218,15 @@ void SimulationMesh::initializeElements() {
2, std::vector<VectorType>(6, VectorType(0, 0, 0)));
element.derivativeT3.resize(
2, std::vector<VectorType>(6, VectorType(0, 0, 0)));
element.derivativeT1_jplus1.resize(6);
element.derivativeT1_j.resize(6);
element.derivativeT1_jplus1.resize(6);
element.derivativeT2_j.resize(6);
element.derivativeT2_jplus1.resize(6);
element.derivativeT3_j.resize(6);
element.derivativeT3_jplus1.resize(6);
element.derivativeR_j.resize(6);
element.derivativeR_jplus1.resize(6);
element.derivativeT1_jplus1.resize(6, VectorType(0, 0, 0));
element.derivativeT1_j.resize(6, VectorType(0, 0, 0));
element.derivativeT1_jplus1.resize(6, VectorType(0, 0, 0));
element.derivativeT2_j.resize(6, VectorType(0, 0, 0));
element.derivativeT2_jplus1.resize(6, VectorType(0, 0, 0));
element.derivativeT3_j.resize(6, VectorType(0, 0, 0));
element.derivativeT3_jplus1.resize(6, VectorType(0, 0, 0));
element.derivativeR_j.resize(6, VectorType(0, 0, 0));
element.derivativeR_jplus1.resize(6, VectorType(0, 0, 0));
}
updateElementalFrames();
}

1
simulationmesh.hpp
View File

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

52
topologyenumerator.cpp
View File

@ -5,7 +5,7 @@
#include <numeric>
#include <unordered_set>
const bool debugIsOn{false};
const bool debugIsOn{true};
const bool savePlyFiles{true};
// size_t binomialCoefficient(size_t n, size_t m) {
@ -456,9 +456,8 @@ void TopologyEnumerator::computeValidPatterns(
.string()
+ ".ply");
}
} else {
continue; // should be uncommented in order to improve performance
}
continue; // should be uncommented in order to improve performance
}
const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold
@ -553,6 +552,50 @@ void TopologyEnumerator::computeValidPatterns(
+ ".ply");
PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(
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)
.append(patternName + "_tiled")
.string()
@ -563,8 +606,7 @@ void TopologyEnumerator::computeValidPatterns(
}
}
if (hasArticulationPoints /*&& !hasFloatingComponents &&
!tiledPatternHasDanglingEdges*/) {
if (hasArticulationPoints && !hasFloatingComponents && !tiledPatternHasDanglingEdges) {
statistics.numberOfPatternsWithArticulationPoints++;
if (debugIsOn) {
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>::EdgeEdge(tile);
for (size_t vi = 0; vi < pattern.vn; vi++) {
tile.vert[vi].C() = vcg::Color4b::Blue;
}
// for (size_t vi = 0; vi < pattern.vn; vi++) {
// tile.vert[vi].C() = vcg::Color4b::Blue;
// }
tile.updateEigenEdgeAndVertices();
return tile;
@ -816,7 +816,7 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
std::vector<ConstPatternGeometry> &patterns,
const std::vector<int> &connectToNeighborsVi,
const VCGPolyMesh &tileInto,
const std::vector<size_t> &perSurfaceFacePatternIndices,
const std::vector<int> &perSurfaceFacePatternIndices,
std::vector<size_t> &tileIntoEdgesToTiledVi,
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::HasFVAdjacency(tileInto));
for (const VCGPolyMesh::FaceType &f : tileInto.face) {
const int facePatternIndex = perSurfaceFacePatternIndices[tileInto.getIndex(f)];
if (facePatternIndex == -1) {
continue;
}
CoordType centerOfFace(0, 0, 0);
for (size_t vi = 0; vi < f.VN(); vi++) {
centerOfFace = centerOfFace + f.cP(vi);
@ -856,7 +860,6 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
for (int &vi : firstInFanConnectToNeighbor_vi) {
vi += pTiledPattern->VN();
}
const size_t facePatternIndex = perSurfaceFacePatternIndices[tileInto.getIndex(f)];
ConstPatternGeometry &pattern = patterns[facePatternIndex];
for (size_t vi = 0; vi < f.VN(); vi++) {
@ -966,6 +969,9 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
tileIntoEdgesToTiledVi.resize(tileInto.EN());
for (int ei = 0; ei < tileInto.EN(); ei++) {
const std::vector<int> interfaceVis = tileIntoEdgeToInterfaceVi[ei];
if (interfaceVis.empty()) {
continue;
}
assert(interfaceVis.size() == 1 || interfaceVis.size() == 2);
assert(
interfaceVis.size() == 1

2
trianglepatterngeometry.hpp
View File

@ -120,7 +120,7 @@ private:
std::vector<ConstPatternGeometry> &patterns,
const std::vector<int> &connectToNeighborsVi,
const VCGPolyMesh &tileInto,
const std::vector<size_t> &perSurfaceFacePatternIndices,
const std::vector<int> &perSurfaceFacePatternIndices,
std::vector<size_t> &tileIntoEdgesToTiledVi,
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;
}
inline size_t computeHashOrdered(const std::vector<size_t> &v)
inline size_t computeHashOrdered(const std::vector<int> &v)
{
std::string elementsString;
for (const auto &el : v) {

35
vcgtrimesh.cpp
View File

@ -8,27 +8,31 @@
bool VCGTriMesh::load(const std::string &filename) {
assert(std::filesystem::exists(filename));
// unsigned int mask = 0;
// mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTCOORD;
// mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTNORMAL;
// mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTCOLOR;
// mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_EDGEINDEX;
// mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_FACEINDEX;
unsigned int mask = 0;
mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTCOORD;
mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTNORMAL;
mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTCOLOR;
mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_EDGEINDEX;
mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_FACEINDEX;
// if (nanoply::NanoPlyWrapper<VCGTriMesh>::LoadModel(
// std::filesystem::absolute(filename).string().c_str(), *this, mask) != 0) {
if (vcg::tri::io::Importer<VCGTriMesh>::Open(*this,
std::filesystem::absolute(filename).string().c_str())
if (nanoply::NanoPlyWrapper<VCGTriMesh>::LoadModel(
std::filesystem::absolute(filename).string().c_str(), *this, mask)
!= 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;
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();
return true;
vcg::tri::UpdateTopology<VCGTriMesh>::AllocateEdge(*this);
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
@ -87,6 +91,7 @@ bool VCGTriMesh::save(const std::string plyFilename)
}
#ifdef POLYSCOPE_DEFINED
polyscope::SurfaceMesh *VCGTriMesh::registerForDrawing(
const std::optional<std::array<double, 3>> &desiredColor, const bool &shouldEnable) const
{

8
vcgtrimesh.hpp
View File

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