diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index ddeadb6..6f3cc15 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -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(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(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> &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 &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::Quaternion 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(); + 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(); + // const VectorType &n_initial = node.initialNormal; + const Eigen::Vector3d referenceF1_initial + = (referenceT1_initial - (nInitial_eigen * (referenceT1_initial.dot(nInitial_eigen)))) + .normalized(); + Eigen::Quaternion q_f1; //nr is with respect to f1 + q_f1.setFromTwoVectors(referenceF1_initial, referenceF1_deformed); + Eigen::Quaternion q_normal; + q_normal.setFromTwoVectors(nInitial_eigen, nDeformed_eigen); + Eigen::Quaternion 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 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 &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 &externalForce : pJob->nodalExternalForces) { - totalExternalForcesNorm += externalForce.second.norm(); - } + reset(); + // double totalExternalForcesNorm = 0; + // for (const std::pair &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_ptrnodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) { - shouldApplyInitialDistortion = true; - } - updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices); - const size_t movingAverageSampleSize = 50; - std::queue 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 residualForcesMovingAverageHistorySample; + std::vector 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(); - const Eigen::Vector3d deformedNormal + const Eigen::Vector3d nInitial_eigen = node.initialNormal.ToEigenVector(); + const Eigen::Vector3d nDeformed_eigen = pMesh->vert[vi].cN().ToEigenVector(); - Eigen::Quaternion q; - q.setFromTwoVectors(initialNormal, deformedNormal); - Eigen::Quaternion q_nr; - q_nr = Eigen::AngleAxis(pMesh->nodes[vi].nR, deformedNormal); + Eigen::Quaternion q_normal; + q_normal.setFromTwoVectors(nInitial_eigen, nDeformed_eigen); + Eigen::Quaternion q_nr_nDeformed; + q_nr_nDeformed = Eigen::AngleAxis(pMesh->nodes[vi].nR, nDeformed_eigen); + Eigen::Quaternion q_nr_nInit; + q_nr_nInit = Eigen::AngleAxis(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 q_f1; //nr is with respect to f1 - q_f1.setFromTwoVectors(referenceF1_initial.ToEigenVector(), - referenceF1_deformed.ToEigenVector()); - Eigen::Quaternion q_t1; - q_t1.setFromTwoVectors(referenceT1_initial.ToEigenVector(), - referenceT1_deformed.ToEigenVector()); + // const VectorType &n_initial = node.initialNormal; + const VectorType referenceF1_initial = (referenceT1_initial + - (node.initialNormal + * (referenceT1_initial * node.initialNormal))) + .Normalize(); + Eigen::Quaternion q_f1_nInit; //nr is with respect to f1 + q_f1_nInit.setFromTwoVectors(referenceF1_initial.ToEigenVector(), + referenceF1_deformed.ToEigenVector()); - results.rotationalDisplacementQuaternion[vi] = q * q_f1 * q_nr; + Eigen::Quaternion 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(), + referenceF1_deformed_def.ToEigenVector()); + 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 q_t1; + // q_t1.setFromTwoVectors(referenceT1_initial.ToEigenVector(), + // referenceT1_deformed.ToEigenVector()); + + 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++; } } diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 4e790fc..392fd95 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -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 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 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 plotYValues; + size_t numOfDampings{0}; - const std::string meshPolyscopeLabel{"Simulation mesh"}; - std::unique_ptr pMesh; - std::unordered_map> - constrainedVertices; - SimulationHistory history; - // Eigen::Tensor theta3Derivatives; - // std::unordered_map theta3Derivatives; - bool shouldApplyInitialDistortion = false; - std::unordered_set rigidSupports; + const std::string meshPolyscopeLabel{"Simulation mesh"}; + std::unique_ptr pMesh; + std::unordered_map> constrainedVertices; + SimulationHistory history; + // Eigen::Tensor theta3Derivatives; + // std::unordered_map theta3Derivatives; + bool shouldApplyInitialDistortion = false; + std::unordered_set rigidSupports; - void reset(); - void updateNodalInternalForces( - const std::unordered_map> - &fixedVertices); - void updateNodalExternalForces( - const std::unordered_map &nodalForces, - const std::unordered_map> - &fixedVertices); - void updateResidualForces(); - void updateRotationalDisplacements(); - void updateElementalLengths(); + void reset(); + void updateNodalInternalForces( + const std::unordered_map> &fixedVertices); + void updateNodalExternalForces( + const std::unordered_map &nodalForces, + const std::unordered_map> &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 &nodalForcedDisplacements); + void applyForcedDisplacements( + const std::unordered_map &nodalForcedDisplacements); - Vector6d computeElementTorsionalForce( - const Element &element, const Vector6d &displacementDifference, - const std::unordered_set &constrainedDof); + Vector6d computeElementTorsionalForce(const Element &element, + const Vector6d &displacementDifference, + const std::unordered_set &constrainedDof); - // BeamFormFinder::Vector6d computeElementFirstBendingForce( - // const Element &element, const Vector6d &displacementDifference, - // const std::unordered_set &constrainedDof); + // BeamFormFinder::Vector6d computeElementFirstBendingForce( + // const Element &element, const Vector6d &displacementDifference, + // const std::unordered_set &constrainedDof); - // BeamFormFinder::Vector6d computeElementSecondBendingForce( - // const Element &element, const Vector6d &displacementDifference, - // const std::unordered_set &constrainedDof); + // BeamFormFinder::Vector6d computeElementSecondBendingForce( + // const Element &element, const Vector6d &displacementDifference, + // const std::unordered_set &constrainedDof); - void updateKineticEnergy(); + void updateKineticEnergy(); - void resetVelocities(); + void resetVelocities(); - SimulationResults computeResults(const std::shared_ptr &pJob); + SimulationResults computeResults(const std::shared_ptr &pJob); - void updateNodePosition( - VertexType &v, - const std::unordered_map> - &fixedVertices); + void updateNodePosition( + VertexType &v, + const std::unordered_map> &fixedVertices); - void applyDisplacements( - const std::unordered_map> - &fixedVertices); + void applyDisplacements( + const std::unordered_map> &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 &pJob); + + void updateNodeNr(VertexType &v); + public: DRMSimulationModel(); - SimulationResults - executeSimulation(const std::shared_ptr &pJob, - const Settings &settings = Settings()); + SimulationResults executeSimulation(const std::shared_ptr &pJob, + const Settings &settings = Settings(), + const SimulationResults &solutionGuess = SimulationResults()); static void runUnitTests(); }; diff --git a/edgemesh.cpp b/edgemesh.cpp index a04813c..2fac507 100755 --- a/edgemesh.cpp +++ b/edgemesh.cpp @@ -220,6 +220,7 @@ bool VCGEdgeMesh::loadUsingNanoply(const std::string &plyFilename) { unsigned int mask = 0; mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; + mask |= nanoply::NanoPlyWrapper::IO_VERTCOLOR; mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; if (nanoply::NanoPlyWrapper::LoadModel(plyFilename.c_str(), *this, mask) != 0) { diff --git a/polymesh.hpp b/polymesh.hpp index 931fca6..ea11c57 100755 --- a/polymesh.hpp +++ b/polymesh.hpp @@ -110,7 +110,7 @@ public: return faces; } #ifdef POLYSCOPE_DEFINED - std::shared_ptr registerForDrawing( + polyscope::SurfaceMesh *registerForDrawing( const std::optional> &desiredColor = std::nullopt, const bool &shouldEnable = true) { @@ -118,8 +118,9 @@ public: auto faces = getFaces(); PolyscopeInterface::init(); - std::shared_ptr 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); diff --git a/reducedmodeloptimizer_structs.hpp b/reducedmodeloptimizer_structs.hpp index c05ab8f..08b92d1 100644 --- a/reducedmodeloptimizer_structs.hpp +++ b/reducedmodeloptimizer_structs.hpp @@ -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(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(job)); + job.load(filepath.string()); + fullPatternSimulationJobs.push_back(std::make_shared(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(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(); diff --git a/simulation_structs.hpp b/simulation_structs.hpp index 895a5dd..e1d3175 100755 --- a/simulation_structs.hpp +++ b/simulation_structs.hpp @@ -35,7 +35,7 @@ struct SimulationHistory { size_t numberOfSteps{0}; std::string label; - std::vector residualForcesLog; + std::vector logResidualForces; std::vector kineticEnergy; std::vector potentialEnergies; std::vector 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 &thisToOtherViMap, SimulationJob &simulationJobMapped) const { @@ -418,14 +424,18 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m ->setEnabled(shouldEnable); } // per node external moments + bool hasExternalMoments = false; std::vector> 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); @@ -436,8 +446,12 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m 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 }; @@ -446,9 +460,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 job; SimulationHistory history; + std::vector debug_drmDisplacements; + std::vector> debug_q_f1; //per vertex + std::vector> debug_q_normal; //per vertex + std::vector> debug_q_nr; //per vertex std::vector displacements; std::vector> rotationalDisplacementQuaternion; //per vertex double internalPotentialEnergy{0}; diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp index 422a4a8..ce317b9 100755 --- a/simulationhistoryplotter.hpp +++ b/simulationhistoryplotter.hpp @@ -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 &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 &XvaluesToPlot, - const std::vector &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 &x, + const std::vector &y, + const std::vector &markerSizes, + const std::vector &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 &YvaluesToPlot, - const std::string &saveTo = {}) + const std::string &saveTo = {}, + const std::vector &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 colors(x.size(), 0.5); + std::vector 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()); diff --git a/simulationmesh.cpp b/simulationmesh.cpp index c2cf67d..997364d 100755 --- a/simulationmesh.cpp +++ b/simulationmesh.cpp @@ -40,20 +40,46 @@ SimulationMesh::SimulationMesh(PatternGeometry &pattern) { initializeElements(); } -SimulationMesh::SimulationMesh(SimulationMesh &mesh) { - vcg::tri::MeshAssert::VertexNormalNormalized(mesh); - VCGEdgeMesh::copy(mesh); +SimulationMesh::SimulationMesh(SimulationMesh &m) +{ + vcg::tri::MeshAssert::VertexNormalNormalized(m); + VCGEdgeMesh::copy(m); - elements = vcg::tri::Allocator::GetPerEdgeAttribute( - *this, std::string("Elements")); - nodes = vcg::tri::Allocator::GetPerVertexAttribute( - *this, std::string("Nodes")); - initializeNodes(); + elements = vcg::tri::Allocator::GetPerEdgeAttribute(*this, + std::string( + "Elements")); + nodes = vcg::tri::Allocator::GetPerVertexAttribute(*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::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::VertexEdge(*this); + + elements = vcg::tri::Allocator::GetPerEdgeAttribute(*this, + std::string( + "Elements")); + nodes = vcg::tri::Allocator::GetPerVertexAttribute(*this, + std::string("Nodes")); + initializeNodes(); + initializeElements(); + reset(); } void SimulationMesh::computeElementalProperties() { @@ -192,15 +218,15 @@ void SimulationMesh::initializeElements() { 2, std::vector(6, VectorType(0, 0, 0))); element.derivativeT3.resize( 2, std::vector(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(); } diff --git a/simulationmesh.hpp b/simulationmesh.hpp index 5dea7af..ca3dc38 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -28,6 +28,7 @@ public: SimulationMesh(PatternGeometry &pattern); SimulationMesh(ConstVCGEdgeMesh &edgeMesh); SimulationMesh(SimulationMesh &elementalMesh); + SimulationMesh(VCGTriMesh &triMesh); void updateElementalLengths(); void updateIncidentElements(); diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp index 0178863..e8c4f99 100755 --- a/topologyenumerator.cpp +++ b/topologyenumerator.cpp @@ -5,7 +5,7 @@ #include #include -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> eCC; + vcg::tri::Clean::edgeMeshConnectedComponents(tiledPatternGeometry, + eCC); + vcg::tri::UpdateFlags::EdgeClear(tiledPatternGeometry); + const size_t numberOfCC_edgeBased = eCC.size(); + std::sort(eCC.begin(), + eCC.end(), + [](const std::pair &a, + const std::pair &b) { + return a.first > b.first; + }); + + PatternGeometry::EdgePointer &ep = eCC[0].second; + size_t colorsRegistered = 0; + std::stack stack; + stack.push(ep); + while (!stack.empty()) { + EdgePointer ep = stack.top(); + stack.pop(); + + for (int i = 0; i < 2; ++i) { + vcg::edge::VEIterator 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) { diff --git a/trianglepatterngeometry.cpp b/trianglepatterngeometry.cpp index e1cbe47..921a768 100755 --- a/trianglepatterngeometry.cpp +++ b/trianglepatterngeometry.cpp @@ -112,9 +112,9 @@ PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern) { vcg::tri::UpdateTopology::VertexEdge(tile); vcg::tri::UpdateTopology::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::tilePattern( std::vector &patterns, const std::vector &connectToNeighborsVi, const VCGPolyMesh &tileInto, - const std::vector &perSurfaceFacePatternIndices, + const std::vector &perSurfaceFacePatternIndices, std::vector &tileIntoEdgesToTiledVi, std::vector> &perPatternIndexToTiledPatternEdgeIndex) { @@ -843,6 +843,10 @@ std::shared_ptr 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::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::tilePattern( tileIntoEdgesToTiledVi.resize(tileInto.EN()); for (int ei = 0; ei < tileInto.EN(); ei++) { const std::vector interfaceVis = tileIntoEdgeToInterfaceVi[ei]; + if (interfaceVis.empty()) { + continue; + } assert(interfaceVis.size() == 1 || interfaceVis.size() == 2); assert( interfaceVis.size() == 1 diff --git a/trianglepatterngeometry.hpp b/trianglepatterngeometry.hpp index 5bbc857..8364830 100755 --- a/trianglepatterngeometry.hpp +++ b/trianglepatterngeometry.hpp @@ -120,7 +120,7 @@ private: std::vector &patterns, const std::vector &connectToNeighborsVi, const VCGPolyMesh &tileInto, - const std::vector &perSurfaceFacePatternIndices, + const std::vector &perSurfaceFacePatternIndices, std::vector &tileIntoEdgesToTiledVi, std::vector> &perPatternIndexTiledPatternEdgeIndex); }; diff --git a/utilities.hpp b/utilities.hpp index d93dac4..cf13769 100755 --- a/utilities.hpp +++ b/utilities.hpp @@ -315,7 +315,7 @@ size_t computeHashUnordered(const std::vector &v) return hash; } -inline size_t computeHashOrdered(const std::vector &v) +inline size_t computeHashOrdered(const std::vector &v) { std::string elementsString; for (const auto &el : v) { diff --git a/vcgtrimesh.cpp b/vcgtrimesh.cpp index c04664a..2f630b7 100755 --- a/vcgtrimesh.cpp +++ b/vcgtrimesh.cpp @@ -8,27 +8,31 @@ bool VCGTriMesh::load(const std::string &filename) { assert(std::filesystem::exists(filename)); - // unsigned int mask = 0; - // mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; - // mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; - // mask |= nanoply::NanoPlyWrapper::IO_VERTCOLOR; - // mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; - // mask |= nanoply::NanoPlyWrapper::IO_FACEINDEX; + unsigned int mask = 0; + mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; + mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; + mask |= nanoply::NanoPlyWrapper::IO_VERTCOLOR; + mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; + mask |= nanoply::NanoPlyWrapper::IO_FACEINDEX; - // if (nanoply::NanoPlyWrapper::LoadModel( - // std::filesystem::absolute(filename).string().c_str(), *this, mask) != 0) { - if (vcg::tri::io::Importer::Open(*this, - std::filesystem::absolute(filename).string().c_str()) + if (nanoply::NanoPlyWrapper::LoadModel( + std::filesystem::absolute(filename).string().c_str(), *this, mask) != 0) { + // if (vcg::tri::io::Importer::Open(*this, + // std::filesystem::absolute(filename).string().c_str()) + // != 0) { std::cout << "Could not load tri mesh" << std::endl; return false; } - vcg::tri::UpdateTopology::FaceFace(*this); - vcg::tri::UpdateTopology::VertexFace(*this); - vcg::tri::UpdateNormal::PerVertexNormalized(*this); - label = std::filesystem::path(filename).stem().string(); - return true; + vcg::tri::UpdateTopology::AllocateEdge(*this); + vcg::tri::UpdateTopology::FaceFace(*this); + vcg::tri::UpdateTopology::VertexFace(*this); + vcg::tri::UpdateTopology::VertexEdge(*this); + vcg::tri::UpdateNormal::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> &desiredColor, const bool &shouldEnable) const { diff --git a/vcgtrimesh.hpp b/vcgtrimesh.hpp index 60d8f9e..23ac9f1 100755 --- a/vcgtrimesh.hpp +++ b/vcgtrimesh.hpp @@ -21,14 +21,15 @@ class VCGTriMeshVertex : public vcg::Vertex + vcg::vertex::Qualityd, + vcg::vertex::VEAdj> {}; class VCGTriMeshFace : public vcg::Face {}; -class VCGTriMeshEdge - : public vcg::Edge {}; +class VCGTriMeshEdge : public vcg::Edge +{}; class VCGTriMesh : public vcg::tri::TriMesh, std::vector, @@ -52,6 +53,7 @@ public: const bool &shouldEnable = true) const; #endif Eigen::MatrixX2i getEdges() const; + void updateEigenEdgeAndVertices(); }; #endif // VCGTRIMESH_HPP