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

Refactoring

This commit is contained in:
iasonmanolas 2021-11-15 11:08:39 +02:00
parent 21d79f07aa
commit f8a98deb17
21 changed files with 1191 additions and 654 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
CMakeLists.txt
View File

@ -16,6 +16,32 @@ if(NOT EXISTS ${EXTERNAL_DEPS_DIR})
endif()
##Create directory for the external libraries
file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR})
#message(${POLYSCOPE_ALREADY_COMPILED})
if(${USE_POLYSCOPE})
download_project(PROJ POLYSCOPE
GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git
GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
if(NOT EXISTS ${POLYSCOPE_BINARY_DIR})
add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR})
add_compile_definitions(POLYSCOPE_DEFINED)
endif()
endif()
#dlib
set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin)
file(MAKE_DIRECTORY ${DLIB_BIN_DIR})
download_project(PROJ DLIB
GIT_REPOSITORY https://github.com/davisking/dlib.git
GIT_TAG master
BINARY_DIR ${DLIB_BIN_DIR}
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR})
##vcglib devel branch
download_project(PROJ vcglib_devel
@ -55,7 +81,6 @@ option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastruct
add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR})
link_directories(${TBB_BINARY_DIR})
###Eigen 3 NOTE: Eigen is required on the system the code is ran
find_package(Eigen3 3.3 REQUIRED)
@ -75,9 +100,11 @@ target_include_directories(${PROJECT_NAME}
PUBLIC ${MySourcesFiles}
)
if(${USE_POLYSCOPE})
target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot polyscope glad ThreedBeamFEA tbb pthread)
if(${MYSOURCES_STATIC_LINK})
target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread)
elseif(${USE_POLYSCOPE})
target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib glad ThreedBeamFEA tbb pthread polyscope)
else()
target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread)
target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib glad ThreedBeamFEA tbb pthread)
endif()
target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)

3
csvfile.hpp
View File

@ -36,8 +36,7 @@ public:
std::ofstream outfile(filename);
outfile.close();
}
overwrite ? fs_.open(filename, std::ios::trunc)
: fs_.open(filename, std::ios::app);
overwrite ? fs_.open(filename, std::ios::trunc) : fs_.open(filename, std::ios::app);
}
}

435
drmsimulationmodel.cpp
View File

@ -218,22 +218,58 @@ void DRMSimulationModel::runUnitTests()
// polyscope::show();
}
void DRMSimulationModel::reset()
void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings)
{
mSettings = settings;
mCurrentSimulationStep = 0;
history.clear();
history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
constrainedVertices.clear();
pMesh.reset();
pMesh = std::make_unique<SimulationMesh>(*pJob->pMesh);
vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
plotYValues.clear();
plotHandle.reset();
checkedForMaximumMoment = false;
externalMomentsNorm = 0;
Dt = mSettings.Dtini;
Dt = settings.Dtini;
numOfDampings = 0;
shouldTemporarilyDampForces = false;
externalLoadStep = 1;
isVertexConstrained.clear();
minTotalResidualForcesNorm = std::numeric_limits<double>::max();
constrainedVertices = pJob->constrainedVertices;
if (!pJob->nodalForcedDisplacements.empty()) {
for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair : pJob->nodalForcedDisplacements) {
const VertexIndex vi = viDisplPair.first;
constrainedVertices[vi].insert({0, 1, 2});
}
}
computeRigidSupports();
isVertexConstrained.resize(pMesh->VN(), false);
for (auto fixedVertex : constrainedVertices) {
isVertexConstrained[fixedVertex.first] = true;
}
#ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) {
PolyscopeInterface::init();
polyscope::registerCurveNetwork(meshPolyscopeLabel,
pMesh->getEigenVertices(),
pMesh->getEigenEdges());
polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel,
pMesh->getEigenVertices(),
pMesh->getEigenEdges())
->setRadius(0.002);
// registerWorldAxes();
}
#endif
if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) {
shouldApplyInitialDistortion = true;
}
updateElementalFrames();
}
VectorType DRMSimulationModel::computeDisplacementDifferenceDerivative(
@ -755,18 +791,18 @@ double DRMSimulationModel::computeTotalInternalPotentialEnergy()
const double axialPE = pow(e_k, 2) * element.material.youngsModulus * element.A
/ (2 * element.initialLength);
const double theta1Diff = theta1_jplus1 - theta1_j;
const double torsionalPE = element.G * element.J * pow(theta1Diff, 2)
const double torsionalPE = element.G * element.inertia.J * pow(theta1Diff, 2)
/ (2 * element.initialLength);
const double firstBendingPE_inBracketsTerm = 4 * pow(theta2_j, 2)
+ 4 * theta2_j * theta2_jplus1
+ 4 * pow(theta2_jplus1, 2);
const double firstBendingPE = firstBendingPE_inBracketsTerm * element.material.youngsModulus
* element.I2 / (2 * element.initialLength);
* element.inertia.I2 / (2 * element.initialLength);
const double secondBendingPE_inBracketsTerm = 4 * pow(theta3_j, 2)
+ 4 * theta3_j * theta3_jplus1
+ 4 * pow(theta3_jplus1, 2);
const double secondBendingPE = secondBendingPE_inBracketsTerm
* element.material.youngsModulus * element.I3
* element.material.youngsModulus * element.inertia.I3
/ (2 * element.initialLength);
totalInternalPotentialEnergy += axialPE + torsionalPE + firstBendingPE + secondBendingPE;
@ -1220,11 +1256,12 @@ void DRMSimulationModel::updateNodalMasses()
const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length;
translationalSumSk += SkTranslational;
const double lengthToThe3 = std::pow(elem.length, 3);
const long double SkRotational_I2 = elem.material.youngsModulus * elem.I2
const long double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2
/ lengthToThe3;
const long double SkRotational_I3 = elem.material.youngsModulus * elem.I3
const long double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3
/ lengthToThe3;
const long double SkRotational_J = elem.material.youngsModulus * elem.J / lengthToThe3;
const long double SkRotational_J = elem.material.youngsModulus * elem.inertia.J
/ lengthToThe3;
rotationalSumSk_I2 += SkRotational_I2;
rotationalSumSk_I3 += SkRotational_I3;
rotationalSumSk_J += SkRotational_J;
@ -1586,11 +1623,11 @@ void DRMSimulationModel::updatePositionsOnTheFly(
const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length;
translationalSumSk += SkTranslational;
const double lengthToThe3 = std::pow(elem.length, 3);
const double SkRotational_I2 = elem.material.youngsModulus * elem.I2
const double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const double SkRotational_I3 = elem.material.youngsModulus * elem.I3
const double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const double SkRotational_J = elem.material.youngsModulus * elem.J
const double SkRotational_J = elem.material.youngsModulus * elem.inertia.J
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
rotationalSumSk_I2 += SkRotational_I2;
rotationalSumSk_I3 += SkRotational_I3;
@ -1686,9 +1723,129 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptr<Simul
history.numberOfSteps = mCurrentSimulationStep;
SimulationResults results;
results.converged = true;
results.job = pJob;
results.pJob = pJob;
results.history = history;
results.displacements = displacements;
results.setLabelPrefix("DRM");
if (mSettings.maxDRMIterations.has_value()
&& mCurrentSimulationStep == mSettings.maxDRMIterations && mCurrentSimulationStep != 0) {
if (mSettings.beVerbose) {
std::cout << "Did not reach equilibrium before reaching the maximum number "
"of DRM steps ("
<< mSettings.maxDRMIterations.value() << "). Breaking simulation"
<< std::endl;
}
results.converged = false;
} else if (std::isnan(pMesh->currentTotalKineticEnergy)) {
if (mSettings.beVerbose) {
std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl;
}
results.converged = false;
}
// mesh.printVertexCoordinates(mesh.VN() / 2);
#ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw && !mSettings.debugModeStep.has_value()) {
draw();
}
#endif
if (!std::isnan(pMesh->currentTotalKineticEnergy)) {
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 nInitial_eigen = node.initialNormal
.ToEigenVector<Eigen::Vector3d>();
const Eigen::Vector3d nDeformed_eigen
= pMesh->vert[vi].cN().ToEigenVector<Eigen::Vector3d>();
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 &nDeformed = pMesh->vert[vi].cN();
const VectorType referenceF1_deformed
= (referenceT1_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
- (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>());
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++;
}
}
return results;
}
@ -1844,9 +2001,9 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
for (const EdgeType &e : pMesh->edge) {
const size_t ei = pMesh->getIndex(e);
A[ei] = pMesh->elements[e].A;
J[ei] = pMesh->elements[e].J;
I2[ei] = pMesh->elements[e].I2;
I3[ei] = pMesh->elements[e].I3;
J[ei] = pMesh->elements[e].inertia.J;
I2[ei] = pMesh->elements[e].inertia.I2;
I3[ei] = pMesh->elements[e].inertia.I3;
}
polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A);
@ -2105,93 +2262,9 @@ void DRMSimulationModel::applyKineticDamping(const std::shared_ptr<SimulationJob
++numOfDampings;
}
SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings,
const SimulationResults &solutionGuess)
SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob)
{
assert(pJob->pMesh.operator bool());
auto beginTime = std::chrono::high_resolution_clock::now();
mSettings = settings;
reset();
history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
// if (!pJob->nodalExternalForces.empty()) {
// double externalForcesNorm = 0;
// for (const auto &externalForce : pJob->nodalExternalForces) {
// externalForcesNorm += externalForce.second.norm();
// }
// mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-2;
// }
constrainedVertices = pJob->constrainedVertices;
if (!pJob->nodalForcedDisplacements.empty()) {
for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair : pJob->nodalForcedDisplacements) {
const VertexIndex vi = viDisplPair.first;
constrainedVertices[vi].insert({0, 1, 2});
}
}
// if (!pJob->nodalForcedNormals.empty()) {
// for (std::pair<VertexIndex, Eigen::Vector3d> viNormalPair :
// pJob->nodalForcedDisplacements) {
// assert(viNormalPair3second[2] >= 0);
// }
// }
pMesh.reset();
pMesh = std::make_unique<SimulationMesh>(*pJob->pMesh);
if (mSettings.beVerbose) {
std::cout << "Executing simulation for mesh with:" << pMesh->VN() << " nodes and "
<< pMesh->EN() << " elements." << std::endl;
}
vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
computeRigidSupports();
isVertexConstrained.resize(pMesh->VN(), false);
for (auto fixedVertex : constrainedVertices) {
isVertexConstrained[fixedVertex.first] = true;
}
#ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) {
PolyscopeInterface::init();
polyscope::registerCurveNetwork(meshPolyscopeLabel,
pMesh->getEigenVertices(),
pMesh->getEigenEdges());
polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel,
pMesh->getEigenVertices(),
pMesh->getEigenEdges())
->setRadius(0.002);
// registerWorldAxes();
}
#endif
if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) {
shouldApplyInitialDistortion = true;
}
updateElementalFrames();
if (!solutionGuess.displacements.empty()) {
assert(!mSettings.linearGuessForceScaleFactor.has_value());
applySolutionGuess(solutionGuess, pJob);
shouldTemporarilyDampForces = true;
}
if (mSettings.linearGuessForceScaleFactor.has_value()) {
const double forceScaleFactor = 1;
for (auto &externalForce : pJob->nodalExternalForces) {
externalForce.second = externalForce.second * forceScaleFactor;
}
LinearSimulationModel linearSimulationModel;
SimulationResults simulationResults_fullPatternLinearModel = linearSimulationModel
.executeSimulation(pJob);
// simulationResults_fullPatternLinearModel.save(fullResultsFolderPath);
for (auto &externalForce : pJob->nodalExternalForces) {
externalForce.second = externalForce.second / forceScaleFactor;
}
applySolutionGuess(simulationResults_fullPatternLinearModel, pJob);
shouldTemporarilyDampForces = true;
}
updateNodalMasses();
// std::unordered_map<VertexIndex, Vector6d> nodalExternalForces = pJob->nodalExternalForces;
// double totalExternalForcesNorm = 0;
@ -2206,7 +2279,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
if (!pJob->nodalExternalForces.empty()) {
mSettings.totalResidualForcesNormThreshold
= settings.totalExternalForcesNormPercentageTermination
= mSettings.totalExternalForcesNormPercentageTermination
* pMesh->totalExternalForcesNorm;
} else {
mSettings.totalResidualForcesNormThreshold = 1e-3;
@ -2220,13 +2293,18 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
<< *mSettings.averageResidualForcesCriterionThreshold << std::endl;
}
}
if (mSettings.beVerbose) {
std::cout << "Executing simulation for mesh with:" << pMesh->VN() << " nodes and "
<< pMesh->EN() << " elements." << std::endl;
}
const size_t movingAverageSampleSize = 200;
std::queue<double> residualForcesMovingAverageHistorySample;
std::vector<double> percentageOfConvergence;
// double residualForcesMovingAverageDerivativeMax = 0;
while (!settings.maxDRMIterations.has_value()
|| mCurrentSimulationStep < settings.maxDRMIterations.value()) {
while (!mSettings.maxDRMIterations.has_value()
|| mCurrentSimulationStep < mSettings.maxDRMIterations.value()) {
if ((mSettings.debugModeStep.has_value() && mCurrentSimulationStep == 50000)) {
// std::filesystem::create_directory("./PatternOptimizationNonConv");
// pJob->save("./PatternOptimizationNonConv");
@ -2532,123 +2610,6 @@ mesh->currentTotalPotentialEnergykN*/
results.executionTime = std::chrono::duration_cast<std::chrono::seconds>(endTime - beginTime)
.count();
if (settings.maxDRMIterations.has_value() && mCurrentSimulationStep == settings.maxDRMIterations
&& mCurrentSimulationStep != 0) {
if (mSettings.beVerbose) {
std::cout << "Did not reach equilibrium before reaching the maximum number "
"of DRM steps ("
<< settings.maxDRMIterations.value() << "). Breaking simulation" << std::endl;
}
results.converged = false;
} else if (std::isnan(pMesh->currentTotalKineticEnergy)) {
if (mSettings.beVerbose) {
std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl;
}
results.converged = false;
}
// mesh.printVertexCoordinates(mesh.VN() / 2);
#ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw && !mSettings.debugModeStep.has_value()) {
draw();
}
#endif
if (!std::isnan(pMesh->currentTotalKineticEnergy)) {
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 nInitial_eigen = node.initialNormal
.ToEigenVector<Eigen::Vector3d>();
const Eigen::Vector3d nDeformed_eigen
= pMesh->vert[vi].cN().ToEigenVector<Eigen::Vector3d>();
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 &nDeformed = pMesh->vert[vi].cN();
const VectorType referenceF1_deformed
= (referenceT1_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
- (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>());
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++;
}
}
if (!mSettings.debugModeStep.has_value() && mSettings.shouldCreatePlots) {
SimulationResultsReporter reporter;
reporter.reportResults({results}, "Results", pJob->pMesh->getLabel());
@ -2660,7 +2621,39 @@ mesh->currentTotalPotentialEnergykN*/
polyscope::removeCurveNetwork("Initial_" + meshPolyscopeLabel);
}
#endif
return results;
}
SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings,
const SimulationResults &solutionGuess)
{
reset(pJob, settings);
assert(pJob->pMesh != nullptr);
if (!solutionGuess.displacements.empty()) {
assert(!mSettings.linearGuessForceScaleFactor.has_value());
applySolutionGuess(solutionGuess, pJob);
shouldTemporarilyDampForces = true;
}
if (mSettings.linearGuessForceScaleFactor.has_value()) {
const double forceScaleFactor = mSettings.linearGuessForceScaleFactor.value();
for (auto &externalForce : pJob->nodalExternalForces) {
externalForce.second = externalForce.second * forceScaleFactor;
}
LinearSimulationModel linearSimulationModel;
SimulationResults simulationResults_fullPatternLinearModel = linearSimulationModel
.executeSimulation(pJob);
// simulationResults_fullPatternLinearModel.save(fullResultsFolderPath);
for (auto &externalForce : pJob->nodalExternalForces) {
externalForce.second = externalForce.second / forceScaleFactor;
}
applySolutionGuess(simulationResults_fullPatternLinearModel, pJob);
shouldTemporarilyDampForces = true;
}
SimulationResults results = executeSimulation(pJob); //calling the virtual function
return results;
}

12
drmsimulationmodel.hpp
View File

@ -1,9 +1,9 @@
#ifndef BEAMFORMFINDER_HPP
#define BEAMFORMFINDER_HPP
#include "simulationmesh.hpp"
#include "matplot/matplot.h"
#include "simulation_structs.hpp"
#include "simulationmesh.hpp"
#include "simulationmodel.hpp"
#include <Eigen/Dense>
#include <filesystem>
#include <unordered_set>
@ -20,7 +20,7 @@ struct DifferentiateWithRespectTo {
const DoFType &dofi;
};
class DRMSimulationModel
class DRMSimulationModel : public SimulationModel
{
public:
struct Settings
@ -95,7 +95,7 @@ private:
// std::unordered_map<MyKeyType, double, key_hash> theta3Derivatives;
bool shouldApplyInitialDistortion = false;
void reset();
void reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings);
void updateNodalInternalForces(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void updateNodalExternalForces(
@ -239,10 +239,12 @@ private:
void applyKineticDamping(const std::shared_ptr<SimulationJob> &pJob);
virtual SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob);
public:
DRMSimulationModel();
SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings = Settings(),
const Settings &settings,
const SimulationResults &solutionGuess = SimulationResults());
static void runUnitTests();

79
edgemesh.cpp
View File

@ -1,9 +1,21 @@
#include "edgemesh.hpp"
#include "vcg/simplex/face/topology.h"
#include <wrap/nanoply/include/nanoplyWrapper.hpp>
//#include <wrap/nanoply/include/nanoplyWrapper.hpp>
#include <wrap/io_trimesh/export.h>
#include <wrap/io_trimesh/import.h>
Eigen::MatrixX2i VCGEdgeMesh::getEigenEdges() const { return eigenEdges; }
std::vector<vcg::Point2i> VCGEdgeMesh::getEdges()
{
getEdges(eigenEdges);
std::vector<vcg::Point2i> edges(eigenEdges.rows());
for (int ei = 0; ei < eigenEdges.rows(); ei++) {
edges[ei] = vcg::Point2i(eigenEdges(ei, 0), eigenEdges(ei, 1));
}
return edges;
}
Eigen::MatrixX3d VCGEdgeMesh::getEigenVertices() {
getVertices(eigenVertices);
return eigenVertices;
@ -23,11 +35,12 @@ bool VCGEdgeMesh::save(const string &plyFilename)
}
assert(std::filesystem::path(filename).extension().string() == ".ply");
unsigned int mask = 0;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
mask |= vcg::tri::io::Mask::IOM_EDGEINDEX;
mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
mask |= vcg::tri::io::Mask::IOM_VERTCOLOR;
if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::SaveModel(filename.c_str(), *this, mask, false) != 0) {
// if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::SaveModel(filename.c_str(), *this, mask, false) != 0) {
if (vcg::tri::io::Exporter<VCGEdgeMesh>::Save(*this, filename.c_str(), mask) != 0) {
return false;
}
return true;
@ -192,9 +205,13 @@ bool VCGEdgeMesh::load(const string &plyFilename) {
}
assert(std::filesystem::exists(usedPath));
this->Clear();
if (!loadUsingNanoply(usedPath)) {
std::cerr << "Error: Unable to open " + usedPath << std::endl;
return false;
// if (!loadUsingNanoply(usedPath)) {
// std::cerr << "Error: Unable to open " + usedPath << std::endl;
// return false;
// }
if (!loadUsingDefaultLoader(usedPath)) {
std::cerr << "Error: Unable to open " + usedPath << std::endl;
return false;
}
getEdges(eigenEdges);
@ -212,21 +229,39 @@ bool VCGEdgeMesh::load(const string &plyFilename) {
return true;
}
bool VCGEdgeMesh::loadUsingNanoply(const std::string &plyFilename) {
//bool VCGEdgeMesh::loadUsingNanoply(const std::string &plyFilename) {
this->Clear();
// assert(plyFileHasAllRequiredFields(plyFilename));
// Load the ply file
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) {
return false;
}
return true;
// this->Clear();
// // assert(plyFileHasAllRequiredFields(plyFilename));
// // Load the ply file
// 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) {
// return false;
// }
// return true;
//}
bool VCGEdgeMesh::loadUsingDefaultLoader(const std::string &plyFilename)
{
this->Clear();
// assert(plyFileHasAllRequiredFields(plyFilename));
// Load the ply file
int mask = 0;
mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
mask |= vcg::tri::io::Mask::IOM_VERTCOLOR;
mask |= vcg::tri::io::Mask::IOM_EDGEINDEX;
// if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::LoadModel(plyFilename.c_str(),
// *this, mask) != 0) {
if (vcg::tri::io::Importer<VCGEdgeMesh>::Open(*this, plyFilename.c_str(), mask) != 0) {
return false;
}
return true;
}
// bool VCGEdgeMesh::plyFileHasAllRequiredFields(const std::string

20
edgemesh.hpp
View File

@ -68,7 +68,7 @@ public:
bool plyFileHasAllRequiredFields(const std::string &plyFilename);
bool loadUsingNanoply(const std::string &plyFilename);
// bool loadUsingNanoply(const std::string &plyFilename);
bool load(const std::string &plyFilename) override;
bool save(const std::string &plyFilename = std::string()) override;
@ -78,13 +78,15 @@ public:
void createSpiral(const float &degreesOfArm, const size_t &numberOfSamples);
Eigen::MatrixX2i getEigenEdges() const;
std::vector<vcg::Point2i> getEdges();
Eigen::MatrixX3d getEigenVertices();
Eigen::MatrixX3d getEigenEdgeNormals() const;
void printVertexCoordinates(const size_t &vi) const;
#ifdef POLYSCOPE_DEFINED
polyscope::CurveNetwork *registerForDrawing(
const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
const double &desiredRadius = 0.002,
const double &desiredRadius = 0.001,
const bool &shouldEnable = true);
void unregister() const;
void drawInitialFrames(polyscope::CurveNetwork *polyscopeHandle_initialMesh) const;
@ -93,10 +95,24 @@ public:
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu_vertices,
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<EdgePointer> &pu_edges);
void centerMesh()
{
CoordType centerOfMass(0, 0, 0);
for (const auto &v : vert) {
centerOfMass = centerOfMass + v.cP();
}
centerOfMass = centerOfMass / VN();
for (auto &v : vert) {
v.P() = v.cP() - centerOfMass;
}
}
private:
void GeneratedRegularSquaredPattern(const double angleDeg,
std::vector<std::vector<vcg::Point2d>> &pattern,
const size_t &desiredNumberOfSamples);
bool loadUsingDefaultLoader(const std::string &plyFilename);
};
using VectorType = VCGEdgeMesh::CoordType;

128
hexagonremesher.hpp Normal file
View File

@ -0,0 +1,128 @@
#ifndef HEXAGONREMESHER_HPP
#define HEXAGONREMESHER_HPP
#include "polymesh.hpp"
#include "vcgtrimesh.hpp"
#include <vcg/complex/algorithms/dual_meshing.h>
#include <vcg/complex/algorithms/isotropic_remeshing.h>
#include <vcg/complex/algorithms/point_sampling.h>
#include <vcg/complex/algorithms/polygon_support.h>
#include <vcg/complex/algorithms/polygonal_algorithms.h>
#include <vcg/complex/algorithms/voronoi_remesher.h>
namespace PolygonalRemeshing {
//std::shared_ptr<VCGPolyMesh>
std::shared_ptr<VCGPolyMesh> remeshWithPolygons(VCGTriMesh &surface)
{
// const std::string tileIntoFilePath(
// "/home/iason/Coding/build/PatternTillingReducedModel/RelWithDebInfo/plane.ply");
// const std::string tileIntoFilePath("/home/iason/Coding/build/PatternTillingReducedModel/"
// "RelWithDebInfo/hexagon.ply");
// const std::string tileIntoFilePath("/home/iason/Coding/build/PatternTillingReducedModel/"
// "RelWithDebInfo/instantMeshes_plane_34.ply");
// const std::string tileIntoFilePath("/home/iason/Coding/build/PatternTillingReducedModel/"
// "RelWithDebInfo/instantMeshes_plane.ply");
// assert(std::filesystem::exists(tileIntoFilePath));
// assert(tileInto_polyMesh.load(tileIntoFilePath));
///Load triangle mesh
// VCGTriMesh tileInto_triMesh;
// assert(std::filesystem::exists(tileIntoFilePath));
// tileInto_triMesh.load(tileIntoFilePath);
///Sample the triangle mesh
// std::vector<CoordType> pointVec;
// double radius = 1;
// int sampleNum = 50;
// // // vcg::tri::PoissonSampling<VCGTriMesh>(tileInto_triMesh, pointVec, sampleNum, radius);
// // // vcg::tri::TrivialSampler<VCGTriMesh> ps;
// // // vcg::tri::SurfaceSampling<VCGTriMesh>::VertexUniform(tileInto_triMesh, ps, sampleNum);
// // // pointVec = ps.SampleVec();
// vcg::tri::MontecarloSampling<VCGTriMesh>(tileInto_triMesh, pointVec, sampleNum);
// vcg::tri::VoronoiProcessingParameter vpp;
// vcg::tri::VoronoiProcessing<VCGTriMesh>::PreprocessForVoronoi(tileInto_triMesh, radius, vpp);
// std::vector<VCGTriMesh::VertexType *> seedVec;
// vcg::tri::VoronoiProcessing<VCGTriMesh>::SeedToVertexConversion(tileInto_triMesh,
// pointVec,
// seedVec);
// vcg::tri::EuclideanDistance<VCGTriMesh> df;
// vpp.geodesicRelaxFlag = false;
// const int iterNum = 100;
// int actualIter = vcg::tri::VoronoiProcessing<VCGTriMesh>::VoronoiRelaxing(tileInto_triMesh,
// seedVec,
// iterNum,
// df,
// vpp);
// std::cout << "Voronoi relaxation iterations performed:" << actualIter << std::endl;
// bool saveWasSuccessful
// = 0
// == vcg::tri::io::ExporterOBJ<VCGTriMesh>::Save(tileInto_triMesh,
// "./tileInto_triMesh.obj",
// vcg::tri::io::Mask::IOM_ALL);
// VCGTriMesh delaunay;
// vcg::tri::VoronoiProcessing<VCGTriMesh>::ConvertDelaunayTriangulationToMesh(tileInto_triMesh,
// delaunay,
// seedVec);
// delaunay.setLabel("Delaunay");
// delaunay.save(delaunay.getLabel() + ".ply");
// delaunay.registerForDrawing();
// ////Load surface
VCGPolyMesh tileInto_polyMesh, dual;
// vcg::tri::PolygonSupport<VCGTriMesh, VCGPolyMesh>::ImportFromTriMesh(tileInto_polyMesh,
// delaunay);
vcg::tri::PolygonSupport<VCGTriMesh, VCGPolyMesh>::ImportFromTriMesh(tileInto_polyMesh, surface);
vcg::tri::DualMeshing<VCGPolyMesh>::MakeDual(tileInto_polyMesh, dual, false);
// for(size_t i=0; i< Edges.size(); ++i)
// {
// std::vector<VCGPolyMesh> fpVec;
// std::vector<int> eiVec;
// face::EFStarFF(Edges[i].f,Edges[i].z,fpVec,eiVec);
// for(size_t j=0;j<fpVec.size();++j)
// fpVec[j]->FEp(eiVec[j])=&(m.edge[i]);
// }
dual.setLabel(surface.getLabel() + "_polyDual");
// bool saveWasSuccessful
// = 0
// == vcg::tri::io::ExporterOBJ<VCGPolyMesh>::Save(dual,
// (dual.getLabel() + ".obj").c_str(),
// vcg::tri::io::Mask::IOM_BITPOLYGONAL);
// assert(saveWasSuccessful);
vcg::tri::UpdateNormal<VCGPolyMesh>::PerFaceNormalized(dual);
// vcg::PolygonalAlgorithm<VCGPolyMesh>::SmoothReprojectPCA(dual, 1000, false, 0.5, 0, false, false);
// vcg::tri::io::ExporterOBJ<VCGPolyMesh>::Save(dual,
// "./dual_optimized.obj",
// vcg::tri::io::Mask::IOM_BITPOLYGONAL);
// if (vcg::tri::VoronoiProcessing<VCGTriMesh>::CheckVoronoiTopology(tileInto_triMesh, seedVec)) {
// VCGTriMesh voronoiMesh, voronoiPoly;
// vcg::tri::VoronoiProcessing<VCGTriMesh>::ConvertVoronoiDiagramToMesh(tileInto_triMesh,
// voronoiMesh,
// voronoiPoly,
// seedVec,
// vpp);
// vcg::tri::Allocator<VCGTriMesh>::CompactEveryVector(voronoiMesh);
// voronoiMesh.setLabel("Voro");
// voronoiMesh.save(voronoiMesh.getLabel() + ".ply");
// // voronoiMesh.registerForDrawing();
// vcg::tri::Allocator<VCGTriMesh>::CompactEveryVector(voronoiPoly);
// voronoiPoly.setLabel("Poly");
// voronoiPoly.save(voronoiPoly.getLabel() + ".ply");
// // voronoiPoly.registerForDrawing();
// // polyscope::show();
// }
// dual.registerForDrawing();
// polyscope::show();
// assert(vcg::tri::IsValidPointer(dual, dual.face[0].cFEp(0)));
return std::make_shared<VCGPolyMesh>(dual);
}
} // namespace PolygonalRemeshing
#endif // HEXAGONREMESHER_HPP

76
linearsimulationmodel.hpp
View File

@ -1,55 +1,13 @@
#ifndef LINEARSIMULATIONMODEL_HPP
#define LINEARSIMULATIONMODEL_HPP
//#include "beam.hpp"
#include "simulation_structs.hpp"
#include "simulationmodel.hpp"
#include "threed_beam_fea.h"
#include <filesystem>
#include <vector>
// struct BeamSimulationProperties {
// float crossArea;
// float I2;
// float I3;
// float polarInertia;
// float G;
// // Properties used by fea
// float EA;
// float EIz;
// float EIy;
// float GJ;
// BeamSimulationProperties(const BeamDimensions &dimensions,
// const BeamMaterial &material);
//};
// struct NodalForce {
// int index;
// int dof;
// double magnitude;
//};
// struct SimulationJob {
// Eigen::MatrixX3d nodes;
// Eigen::MatrixX2i elements;
// Eigen::MatrixX3d elementalNormals;
// Eigen::VectorXi fixedNodes;
// std::vector<NodalForce> nodalForces;
// std::vector<BeamDimensions> beamDimensions;
// std::vector<BeamMaterial> beamMaterial;
//};
// struct SimulationResults {
// std::vector<Eigen::VectorXd> edgeForces; ///< Force values per force
// component
// ///< #force components x #edges
// Eigen::MatrixXd
// nodalDisplacements; ///< The displacement of each node #nodes x 3
// SimulationResults(const fea::Summary &feaSummary);
// SimulationResults() {}
//};
class LinearSimulationModel {
class LinearSimulationModel : public SimulationModel
{
public:
LinearSimulationModel(){
@ -67,8 +25,11 @@ public:
(evn0[2] + evn1[2]) / 2};
const Element &element = job->pMesh->elements[edgeIndex];
const double E = element.material.youngsModulus;
fea::Props feaProperties(E * element.A, E * element.I3, E * element.I2,
element.G * element.J, nAverageVector);
fea::Props feaProperties(E * element.A,
E * element.inertia.I3,
E * element.inertia.I2,
element.G * element.inertia.J,
nAverageVector);
const int vi0 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(0));
const int vi1 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(1));
elements[edgeIndex] = fea::Elem(vi0, vi1, feaProperties);
@ -91,8 +52,8 @@ public:
static std::vector<fea::BC>
getFeaFixedNodes(const std::shared_ptr<SimulationJob> &job) {
std::vector<fea::BC> boundaryConditions;
boundaryConditions.reserve(job->constrainedVertices.size() * 6
+ job->nodalForcedDisplacements.size() * 3);
// boundaryConditions.reserve(job->constrainedVertices.size() * 6
// + job->nodalForcedDisplacements.size() * 3);
for (auto fixedVertex : job->constrainedVertices) {
const int vertexIndex = fixedVertex.first;
for (int dofIndex : fixedVertex.second) {
@ -157,6 +118,7 @@ public:
// results.rotationalDisplacementQuaternion[vi] = q_nz * q_ny * q_nx;
// results.rotationalDisplacementQuaternion[vi] = q_nz * q_nx * q_ny;
}
results.setLabelPrefix("Linear");
// // Convert forces
// // Convert to vector of eigen matrices of the form force component-> per
@ -212,19 +174,19 @@ public:
const double elementPotentialEnergy_bendingY_v0 = std::pow(elementForces[4], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.I2);
* element.inertia.I2);
const double elementPotentialEnergy_bendingZ_v0 = std::pow(elementForces[5], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.I3);
* element.inertia.I3);
const double elementPotentialEnergy_bendingY_v1 = std::pow(elementForces[10], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.I2);
* element.inertia.I2);
const double elementPotentialEnergy_bendingZ_v1 = std::pow(elementForces[11], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.I3);
* element.inertia.I3);
const double elementPotentialEnergy_bending = elementPotentialEnergy_bendingY_v0
+ elementPotentialEnergy_bendingZ_v0
+ elementPotentialEnergy_bendingY_v1
@ -232,10 +194,10 @@ public:
//Torsion
const double elementPotentialEnergy_torsion_v0 = std::pow(elementForces[3], 2)
* element.initialLength
/ (4 * element.J * element.G);
/ (4 * element.inertia.J * element.G);
const double elementPotentialEnergy_torsion_v1 = std::pow(elementForces[9], 2)
* element.initialLength
/ (4 * element.J * element.G);
/ (4 * element.inertia.J * element.G);
const double elementPotentialEnergy_torsion = elementPotentialEnergy_torsion_v0
+ elementPotentialEnergy_torsion_v1;
const double elementInternalPotentialEnergy = elementPotentialEnergy_axial
@ -244,7 +206,7 @@ public:
+ elementPotentialEnergy_torsion;
results.internalPotentialEnergy += elementInternalPotentialEnergy;
}
results.job = simulationJob;
results.pJob = simulationJob;
return results;
}
@ -285,7 +247,7 @@ public:
fea::solve(job, fixedVertices, nodalForces, ties, equations, opts);
SimulationResults results = getResults(feaResults, simulationJob);
results.job = simulationJob;
results.pJob = simulationJob;
return results;
}
// SimulationResults getResults() const;

51
polymesh.hpp
View File

@ -3,7 +3,8 @@
#include "mesh.hpp"
#include "vcg/complex/complex.h"
#include <filesystem>
#include <wrap/io_trimesh/import.h>
#include <wrap/io_trimesh/export_obj.h>
#include <wrap/io_trimesh/export_ply.h>
#ifdef POLYSCOPE_DEFINED
#include <polyscope/surface_mesh.h>
@ -54,7 +55,8 @@ class PFace : public vcg::Face<PUsedTypes,
vcg::face::Qualityf // quality
,
vcg::face::Normal3f // normal
>
,
vcg::face::Color4b>
{};
class VCGPolyMesh
@ -110,6 +112,49 @@ public:
return faces;
}
bool saveOBJ(const std::filesystem::path &objFilePath = std::filesystem::path())
{
const std::string extension = ".obj";
std::filesystem::path filePath = objFilePath;
if (filePath.empty()) {
filePath = std::filesystem::current_path().append(getLabel() + extension).string();
} else if (std::filesystem::is_directory(std::filesystem::path(objFilePath))) {
filePath = std::filesystem::path(objFilePath).append(getLabel() + extension).string();
}
assert(std::filesystem::path(filePath).extension().string() == extension);
unsigned int mask = 0;
mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
if (vcg::tri::io::ExporterOBJ<VCGPolyMesh>::Save(*this, filePath.c_str(), mask) != 0) {
return false;
}
return true;
}
bool savePLY(const std::filesystem::path &objFilePath = std::filesystem::path())
{
const std::string extension = ".ply";
std::filesystem::path filePath = objFilePath;
if (filePath.empty()) {
filePath = std::filesystem::current_path().append(getLabel() + extension).string();
} else if (std::filesystem::is_directory(std::filesystem::path(objFilePath))) {
filePath = std::filesystem::path(objFilePath).append(getLabel() + extension).string();
}
assert(std::filesystem::path(filePath).extension().string() == extension);
unsigned int mask = 0;
mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
if (vcg::tri::io::ExporterPLY<VCGPolyMesh>::Save(*this, filePath.c_str(), mask, false) != 0) {
return false;
}
return true;
}
#ifdef POLYSCOPE_DEFINED
polyscope::SurfaceMesh *registerForDrawing(
const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
@ -131,7 +176,7 @@ public:
const glm::vec3 desiredColor_glm(desiredColor.value()[0],
desiredColor.value()[1],
desiredColor.value()[2]);
polyscopeHandle_mesh->setSurfaceColor(desiredColor_glm);
polyscopeHandle_mesh->setSurfaceColor(glm::normalize(desiredColor_glm));
}
return polyscopeHandle_mesh;

254
reducedmodeloptimizer_structs.hpp
View File

@ -25,19 +25,23 @@ inline static std::vector<std::string> baseSimulationScenarioNames{"Axial",
struct Colors
{
inline static std::array<double, 3> fullInitial{0.416666, 0.109804, 0.890196};
inline static std::array<double, 3> fullDeformed{0.583333, 0.890196, 0.109804};
inline static std::array<double, 3> fullInitial{0.094, 0.094, 0.094};
// inline static std::array<double, 3> fullDeformed{0.583333, 0.890196, 0.109804};
inline static std::array<double, 3> fullDeformed{0.094, 0.094, 0.094};
inline static std::array<double, 3> reducedInitial{0.890196, 0.109804, 0.193138};
inline static std::array<double, 3> reducedDeformed{0.109804, 0.890196, 0.806863};
inline static std::array<double, 3> reducedDeformed{0.262, 0.627, 0.910};
};
struct xRange{
std::string label;
struct xRange
{
std::string label;
double min;
double max;
std::string toString() const {
return label + "=[" + std::to_string(min) + "," + std::to_string(max) +
"]";
inline bool operator<(const xRange &other) { return label < other.label; }
std::string toString() const
{
return label + "=[" + std::to_string(min) + "," + std::to_string(max) + "]";
}
void fromString(const std::string &s)
@ -56,7 +60,8 @@ struct Colors
{
return label == xrange.label && min == xrange.min && max == xrange.max;
}
};
};
struct Settings
{
enum NormalizationStrategy { NonNormalized, Epsilon };
@ -65,7 +70,9 @@ struct Colors
int numberOfFunctionCalls{100000};
double solverAccuracy{1e-2};
NormalizationStrategy normalizationStrategy{Epsilon};
double normalizationParameter{0.003};
double translationNormalizationParameter{0.003};
double rotationNormalizationParameter{3};
bool splitGeometryMaterialOptimization{false};
struct ObjectiveWeights
{
double translational{1};
@ -168,6 +175,7 @@ struct Colors
os << "Normalization strategy";
os << "Trans weight";
os << "Rot weight";
os << "Splitted geo from mat opt";
// os << std::endl;
}
@ -182,9 +190,10 @@ struct Colors
os << numberOfFunctionCalls;
os << solverAccuracy;
os << normalizationStrategyStrings[normalizationStrategy] + "_"
+ std::to_string(normalizationParameter);
+ std::to_string(translationNormalizationParameter);
os << objectiveWeights.translational;
os << objectiveWeights.rotational;
os << (splitGeometryMaterialOptimization == true ? "true" : "false");
}
};
@ -193,8 +202,9 @@ struct Colors
return settings1.numberOfFunctionCalls == settings2.numberOfFunctionCalls
&& settings1.xRanges == settings2.xRanges
&& settings1.solverAccuracy == settings2.solverAccuracy
&& settings1.objectiveWeights.translational
== settings2.objectiveWeights.translational;
&& settings1.objectiveWeights.translational == settings2.objectiveWeights.translational
&& settings1.translationNormalizationParameter
== settings2.translationNormalizationParameter;
}
inline void updateMeshWithOptimalVariables(const std::vector<double> &x, SimulationMesh &m)
@ -212,9 +222,9 @@ struct Colors
Element &e = m.elements[ei];
e.setDimensions(CrossSectionType(beamWidth, beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
e.J = J;
e.I2 = I2;
e.I3 = I3;
e.inertia.J = J;
e.inertia.I2 = I2;
e.inertia.I3 = I3;
}
CoordType center_barycentric(1, 0, 0);
@ -257,6 +267,16 @@ struct Colors
bool wasSuccessful{true};
// int numberOfSimulationCrashes{0};
Settings settings;
std::vector<std::pair<std::string, double>> optimalXNameValuePairs;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern
vcg::Triangle3<double> baseTriangle;
std::string notes;
//Data gathered for csv exporting
struct ObjectiveValues
{
double totalRaw;
@ -268,15 +288,18 @@ struct Colors
std::vector<double> perSimulationScenario_total;
} objectiveValue;
std::vector<double> perScenario_fullPatternPotentialEnergy;
// std::vector<std::pair<std::string,double>> optimalXNameValuePairs;
std::vector<std::pair<std::string, double>> optimalXNameValuePairs;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
std::vector<double> objectiveValueHistory;
std::vector<size_t> objectiveValueHistory_iteration;
PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern
vcg::Triangle3<double> baseTriangle;
std::string notes;
struct CSVExportingSettings
{
bool exportPE{false};
bool exportIterationOfMinima{true};
bool exportRawObjectiveValue{false};
CSVExportingSettings() {}
};
const CSVExportingSettings exportSettings;
struct JsonKeys
{
inline static std::string filename{"OptimizationResults.json"};
@ -290,10 +313,14 @@ struct Colors
void save(const string &saveToPath, const bool shouldExportDebugFiles = false)
{
std::filesystem::create_directories(saveToPath);
//clear directory
for (const auto &entry : std::filesystem::directory_iterator(saveToPath))
std::filesystem::remove_all(entry.path());
if (std::filesystem::exists(saveToPath)) {
for (const auto &entry : std::filesystem::directory_iterator(saveToPath)) {
std::error_code ec;
std::filesystem::remove_all(entry.path(), ec);
}
}
std::filesystem::create_directories(saveToPath);
//Save optimal X
nlohmann::json json_optimizationResults;
@ -373,6 +400,12 @@ struct Colors
}
}
//save minima info
std::filesystem::path csvFilepathMinimaInfo = std::filesystem::path(saveToPath)
.append("minimaInfo.csv");
csvFile csv_minimaInfo(csvFilepathMinimaInfo, false);
writeMinimaInfoTo(csv_minimaInfo);
settings.save(saveToPath);
#ifdef POLYSCOPE_DEFINED
@ -460,7 +493,7 @@ struct Colors
if (filepath.extension() == ".json") {
SimulationJob job;
job.load(filepath.string());
applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh);
// applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh);
applyOptimizationResults_elements(*this, job.pMesh);
reducedPatternSimulationJobs.push_back(
std::make_shared<SimulationJob>(job));
@ -511,6 +544,15 @@ struct Colors
vcg::math::ToRad(hexAngle))
* reducedPattern.vert[0].cP();
}
// for (int rotationCounter = 0;
// rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {
// pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
// = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
// vcg::math::ToRad(60.0 * rotationCounter))
// * baseTriangleMovableVertexPosition;
// }
// reducedPattern.registerForDrawing();
// polyscope::show();
// CoordType p0 = reducedPattern.vert[0].P();
// CoordType p1 = reducedPattern.vert[1].P();
// int i = 0;
@ -528,36 +570,53 @@ struct Colors
reducedPattern_optimizationResults.optimalXNameValuePairs.end());
const std::string ALabel = "A";
assert(optimalXVariables.contains(ALabel));
const double A = optimalXVariables.at(ALabel);
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
CrossSectionType elementDimensions(beamWidth, beamHeight);
if (optimalXVariables.contains(ALabel)) {
const double A = optimalXVariables.at(ALabel);
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
CrossSectionType elementDimensions(beamWidth, beamHeight);
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.setDimensions(elementDimensions);
}
}
const double poissonsRatio = 0.3;
const std::string ymLabel = "E";
assert(optimalXVariables.contains(ymLabel));
const double E = optimalXVariables.at(ymLabel);
const ElementMaterial elementMaterial(poissonsRatio, E);
if (optimalXVariables.contains(ymLabel)) {
const double E = optimalXVariables.at(ymLabel);
const ElementMaterial elementMaterial(poissonsRatio, E);
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.setMaterial(elementMaterial);
}
}
const std::string JLabel = "J";
assert(optimalXVariables.contains(JLabel));
const double J = optimalXVariables.at(JLabel);
if (optimalXVariables.contains(JLabel)) {
const double J = optimalXVariables.at(JLabel);
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.inertia.J = J;
}
}
const std::string I2Label = "I2";
assert(optimalXVariables.contains(I2Label));
const double I2 = optimalXVariables.at(I2Label);
if (optimalXVariables.contains(I2Label)) {
const double I2 = optimalXVariables.at(I2Label);
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.inertia.I2 = I2;
}
}
const std::string I3Label = "I3";
assert(optimalXVariables.contains(I3Label));
const double I3 = optimalXVariables.at(I3Label);
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.setDimensions(elementDimensions);
e.setMaterial(elementMaterial);
e.J = J;
e.I2 = I2;
e.I3 = I3;
if (optimalXVariables.contains(I3Label)) {
const double I3 = optimalXVariables.at(I3Label);
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.inertia.I3 = I3;
}
}
pTiledReducedPattern_simulationMesh->reset();
}
@ -579,9 +638,10 @@ struct Colors
= fullPatternSimulationJobs[simulationJobIndex];
pFullPatternSimulationJob->registerForDrawing(
fullPatternSimulationJobs[0]->pMesh->getLabel());
SimulationResults fullModelResults = drmSimulator.executeSimulation(
pFullPatternSimulationJob);
fullModelResults.registerForDrawing(Colors::fullDeformed, true, true);
DRMSimulationModel::Settings drmSettings;
SimulationResults fullModelResults
= drmSimulator.executeSimulation(pFullPatternSimulationJob, drmSettings);
fullModelResults.registerForDrawing(Colors::fullDeformed, true);
// SimulationResults fullModelLinearResults =
// linearSimulator.executeSimulation(pFullPatternSimulationJob);
// fullModelLinearResults.setLabelPrefix("linear");
@ -589,15 +649,20 @@ struct Colors
// Drawing of reduced pattern results
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob
= reducedPatternSimulationJobs[simulationJobIndex];
// SimulationResults reducedModelResults = drmSimulator.executeSimulation(
// pReducedPatternSimulationJob);
// pReducedPatternSimulationJob->pMesh->registerForDrawing();
// polyscope::show();
// SimulationResults reducedModelResults = drmSimulator.executeSimulation(
// pReducedPatternSimulationJob);
// reducedModelResults.registerForDrawing(Colors::reducedDeformed, false);
SimulationResults reducedModelLinearResults = linearSimulator.executeSimulation(
// SimulationResults reducedModelResults
// = drmSimulator.executeSimulation(pReducedPatternSimulationJob,
// DRMSimulationModel::Settings());
SimulationResults reducedModelResults = linearSimulator.executeSimulation(
pReducedPatternSimulationJob);
reducedModelLinearResults.setLabelPrefix("linear");
reducedModelLinearResults.registerForDrawing(Colors::reducedDeformed, false);
reducedModelResults.setLabelPrefix("linear");
reducedModelResults.registerForDrawing(Colors::reducedDeformed, true);
polyscope::options::programName = fullPatternSimulationJobs[0]->pMesh->getLabel();
polyscope::view::resetCameraToHomeView();
// polyscope::view::resetCameraToHomeView();
polyscope::show();
// Save a screensh
const std::string screenshotFilename
@ -606,18 +671,16 @@ struct Colors
+ fullPatternSimulationJobs[0]->pMesh->getLabel() + "_"
+ pFullPatternSimulationJob->getLabel();
polyscope::screenshot(screenshotFilename, false);
pFullPatternSimulationJob->unregister(fullPatternSimulationJobs[0]->pMesh->getLabel());
fullModelResults.unregister();
// reducedModelResults.unregister();
reducedModelLinearResults.unregister();
reducedModelResults.unregister();
// fullModelLinearResults.unregister();
// double error = computeError(
// reducedModelResults.displacements,fullModelResults.displacements,
// );
// std::cout << "Error of simulation scenario "
// <<
// simula simulationScenarioStrings[simulationScenarioIndex]
// << " is "
// << error << std::endl;
// double error = computeError(reducedModelLinearResults.displacements,
// fullModelResults.displacements);
// std::cout << "Error of simulation scenario "
// << simulationScenarioStrings[simulationScenarioIndex] << " is " << error
// << std::endl;
}
}
#endif // POLYSCOPE_DEFINED
@ -638,8 +701,9 @@ struct Colors
// Drawing of full pattern results
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
= fullPatternSimulationJobs[simulationJobIndex];
SimulationResults fullModelResults = simulator_drm.executeSimulation(
pFullPatternSimulationJob);
DRMSimulationModel::Settings drmSettings;
SimulationResults fullModelResults
= simulator_drm.executeSimulation(pFullPatternSimulationJob, drmSettings);
fullModelResults.saveDeformedModel();
// Drawing of reduced pattern results
@ -653,8 +717,13 @@ struct Colors
void writeHeaderTo(csvFile &os)
{
os << "Total raw Obj value";
if (exportSettings.exportRawObjectiveValue) {
os << "Total raw Obj value";
}
os << "Total Obj value";
if (exportSettings.exportIterationOfMinima) {
os << "Iteration of minima";
}
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < fullPatternSimulationJobs.size();
simulationScenarioIndex++) {
@ -664,12 +733,15 @@ struct Colors
os << "Obj value Rot " + simulationScenarioName;
os << "Obj value Total " + simulationScenarioName;
}
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < fullPatternSimulationJobs.size();
simulationScenarioIndex++) {
const std::string simulationScenarioName
= fullPatternSimulationJobs[simulationScenarioIndex]->getLabel();
os << "PE " + simulationScenarioName;
if (exportSettings.exportPE) {
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < fullPatternSimulationJobs.size();
simulationScenarioIndex++) {
const std::string simulationScenarioName
= fullPatternSimulationJobs[simulationScenarioIndex]->getLabel();
os << "PE " + simulationScenarioName;
}
}
for (const auto &nameValuePair : optimalXNameValuePairs) {
os << nameValuePair.first;
@ -679,10 +751,15 @@ struct Colors
// os << "notes";
}
void writeResultsTo(const Settings &settings_optimization, csvFile &os) const
void writeResultsTo(csvFile &os) const
{
os << objectiveValue.totalRaw;
if (exportSettings.exportRawObjectiveValue) {
os << objectiveValue.totalRaw;
}
os << objectiveValue.total;
if (exportSettings.exportIterationOfMinima) {
os << objectiveValueHistory_iteration.back();
}
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < fullPatternSimulationJobs.size();
simulationScenarioIndex++) {
@ -690,10 +767,12 @@ struct Colors
os << objectiveValue.perSimulationScenario_rotational[simulationScenarioIndex];
os << objectiveValue.perSimulationScenario_total[simulationScenarioIndex];
}
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < fullPatternSimulationJobs.size();
simulationScenarioIndex++) {
os << perScenario_fullPatternPotentialEnergy[simulationScenarioIndex];
if (exportSettings.exportPE) {
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < fullPatternSimulationJobs.size();
simulationScenarioIndex++) {
os << perScenario_fullPatternPotentialEnergy[simulationScenarioIndex];
}
}
for (const auto &optimalXNameValuePair : optimalXNameValuePairs) {
os << optimalXNameValuePair.second;
@ -707,7 +786,20 @@ struct Colors
// }
// os << notes;
}
void writeMinimaInfoTo(csvFile &outputCsv)
{
outputCsv << "Iteration";
outputCsv << "Objective value";
for (int objectiveValueIndex = 0; objectiveValueIndex < objectiveValueHistory.size();
objectiveValueIndex++) {
outputCsv.endrow();
outputCsv << objectiveValueHistory_iteration[objectiveValueIndex];
outputCsv << objectiveValueHistory[objectiveValueIndex];
}
}
};
enum SimulationModelTag { DRM, Linear };
} // namespace ReducedPatternOptimization
#endif // REDUCEDMODELOPTIMIZER_STRUCTS_HPP

156
simulation_structs.hpp
View File

@ -23,6 +23,13 @@ void read_binary(const std::string &filename, Matrix &matrix) {
in.read((char *)matrix.data(), rows * cols * sizeof(typename Matrix::Scalar));
in.close();
}
const static IOFormat CSVFormat(StreamPrecision, DontAlignCols, ", ", "\n");
template<class Matrix>
void writeToCSV(const std::string &filename, Matrix &matrix)
{
ofstream file(filename.c_str());
file << matrix.format(CSVFormat);
}
} // namespace Eigen
#include "simulationmesh.hpp"
@ -195,7 +202,6 @@ public:
bool load(const std::string &jsonFilename, const bool &shouldLoadMesh = true)
{
std::cout << jsonFilename << std::endl;
label = "empty_job";
constrainedVertices.clear();
nodalExternalForces.clear();
@ -406,11 +412,11 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
}
});
if (!nodeColors.empty()) {
polyscope::getCurveNetwork(meshLabel)
->addNodeColorQuantity("Boundary conditions_" + label, nodeColors)
->setEnabled(shouldEnable);
}
// if (!nodeColors.empty()) {
// polyscope::getCurveNetwork(meshLabel)
// ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors)
// ->setEnabled(shouldEnable);
// }
// per node external forces
std::vector<std::array<double, 3>> externalForces(pMesh->VN());
@ -421,9 +427,15 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
}
if (!externalForces.empty()) {
polyscope::getCurveNetwork(meshLabel)
->addNodeVectorQuantity("External force_" + label, externalForces)
->setEnabled(shouldEnable);
polyscope::CurveNetworkNodeVectorQuantity *externalForcesVectors
= polyscope::getCurveNetwork(meshLabel)->addNodeVectorQuantity("External force_"
+ label,
externalForces);
const std::array<double, 3> color_loads{1.0, 0, 0};
externalForcesVectors->setVectorColor(
glm::vec3(color_loads[0], color_loads[1], color_loads[2]));
externalForcesVectors->setEnabled(shouldEnable);
}
// per node external moments
bool hasExternalMoments = false;
@ -454,6 +466,19 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
if (!constrainedVertices.empty()) {
polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label);
}
// per node external moments
bool hasExternalMoments = false;
for (const auto &forcePair : nodalExternalForces) {
const Vector6d &load = forcePair.second;
if (load.getRotation().norm() != 0) {
hasExternalMoments = true;
break;
}
}
if (hasExternalMoments) {
polyscope::getCurveNetwork(meshLabel)->removeQuantity("External moment_" + label);
}
}
#endif // POLYSCOPE_DEFINED
};
@ -463,7 +488,7 @@ struct SimulationResults
* vector contains the same info using euler angles
*/
bool converged{false};
std::shared_ptr<SimulationJob> job;
std::shared_ptr<SimulationJob> pJob;
SimulationHistory history;
std::vector<Vector6d> debug_drmDisplacements;
std::vector<Eigen::Quaternion<double>> debug_q_f1; //per vertex
@ -475,7 +500,7 @@ struct SimulationResults
double executionTime{0};
std::string labelPrefix{"deformed"};
inline static char deliminator{' '};
SimulationResults() { job = std::make_shared<SimulationJob>(); }
SimulationResults() { pJob = std::make_shared<SimulationJob>(); }
std::vector<VectorType> getTranslationalDisplacements() const
{
@ -491,13 +516,13 @@ struct SimulationResults
void setLabelPrefix(const std::string &lp) { labelPrefix += deliminator + lp; }
std::string getLabel() const
{
return labelPrefix + deliminator + job->pMesh->getLabel() + deliminator + job->getLabel();
return labelPrefix + deliminator + pJob->pMesh->getLabel() + deliminator + pJob->getLabel();
}
bool saveDeformedModel(const std::string &outputFolder = std::string())
{
VCGEdgeMesh m;
vcg::tri::Append<VCGEdgeMesh, SimulationMesh>::MeshCopy(m, *job->pMesh);
vcg::tri::Append<VCGEdgeMesh, SimulationMesh>::MeshCopy(m, *pJob->pMesh);
for (int vi = 0; vi < m.VN(); vi++) {
m.vert[vi].P() = m.vert[vi].P()
@ -513,10 +538,11 @@ struct SimulationResults
const std::filesystem::path outputFolderPath = outputFolder.empty()
? std::filesystem::current_path()
: std::filesystem::path(outputFolder);
std::cout << "Saving results to:" << outputFolderPath << std::endl;
std::filesystem::path simulationJobOutputFolderPath
= std::filesystem::path(outputFolderPath).append("SimulationJob");
std::filesystem::create_directories(simulationJobOutputFolderPath);
job->save(simulationJobOutputFolderPath.string());
pJob->save(simulationJobOutputFolderPath.string());
const std::string filename(getLabel() + "_displacements.eigenBin");
Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements);
@ -543,7 +569,7 @@ struct SimulationResults
void load(const std::filesystem::path &loadFromPath, const std::filesystem::path &loadJobFrom)
{
//load job
job->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json").string());
pJob->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json").string());
//Use the first .eigenBin file for loading the displacements
for (auto const &entry : std::filesystem::recursive_directory_iterator(loadFromPath)) {
if (filesystem::is_regular_file(entry) && entry.path().extension() == ".eigenBin") {
@ -554,8 +580,8 @@ struct SimulationResults
}
}
rotationalDisplacementQuaternion.resize(job->pMesh->VN());
for (int vi = 0; vi < job->pMesh->VN(); vi++) {
rotationalDisplacementQuaternion.resize(pJob->pMesh->VN());
for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
rotationalDisplacementQuaternion[vi] = Eigen::AngleAxisd(displacements[vi][3],
Eigen::Vector3d::UnitX())
* Eigen::AngleAxisd(displacements[vi][4],
@ -573,15 +599,18 @@ struct SimulationResults
std::cerr << "Nothing to remove." << std::endl;
return;
}
job->unregister(getLabel());
pJob->unregister(getLabel());
polyscope::removeCurveNetwork(getLabel());
}
void registerForDrawing(const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
const bool &shouldEnable = true,
const bool &shouldShowFrames = false) const
polyscope::CurveNetwork *registerForDrawing(
const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
const bool &shouldEnable = true,
const double &desiredRadius = 0.001,
// const double &desiredRadius = 0.0001,
const bool &shouldShowFrames = false) const
{
PolyscopeInterface::init();
const std::shared_ptr<SimulationMesh> &mesh = job->pMesh;
const std::shared_ptr<SimulationMesh> &mesh = pJob->pMesh;
polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh;
if (!polyscope::hasStructure(getLabel())) {
polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(),
@ -589,7 +618,7 @@ struct SimulationResults
mesh->getEigenEdges());
}
polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable);
polyscopeHandle_deformedEdmeMesh->setRadius(0.002, true);
polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, true);
if (desiredColor.has_value()) {
const glm::vec3 desiredColor_glm(desiredColor.value()[0],
desiredColor.value()[1],
@ -600,42 +629,94 @@ struct SimulationResults
Eigen::MatrixX3d framesX(mesh->VN(), 3);
Eigen::MatrixX3d framesY(mesh->VN(), 3);
Eigen::MatrixX3d framesZ(mesh->VN(), 3);
Eigen::MatrixX3d framesX_initial(mesh->VN(), 3);
Eigen::MatrixX3d framesY_initial(mesh->VN(), 3);
Eigen::MatrixX3d framesZ_initial(mesh->VN(), 3);
// std::unordered_set<int> interfaceNodes{1, 3, 5, 7, 9, 11};
std::unordered_set<int> interfaceNodes{3, 7, 11, 15, 19, 23};
for (VertexIndex vi = 0; vi < mesh->VN(); vi++) {
const Vector6d &nodalDisplacement = displacements[vi];
nodalDisplacements.row(vi) = Eigen::Vector3d(nodalDisplacement[0],
nodalDisplacement[1],
nodalDisplacement[2]);
auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1);
auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0);
auto deformedFrameX = rotationalDisplacementQuaternion[vi] * 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]);
if (interfaceNodes.contains(vi)) {
auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1);
auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0);
auto deformedFrameX = rotationalDisplacementQuaternion[vi] * 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]);
framesX_initial.row(vi) = Eigen::Vector3d(1, 0, 0);
framesY_initial.row(vi) = Eigen::Vector3d(0, 1, 0);
framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 1);
} else {
framesX.row(vi) = Eigen::Vector3d(0, 0, 0);
framesY.row(vi) = Eigen::Vector3d(0, 0, 0);
framesZ.row(vi) = Eigen::Vector3d(0, 0, 0);
framesX_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
}
}
polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices()
+ nodalDisplacements);
const double frameRadius_default = 0.035;
const double frameLength_default = 0.035;
const bool shouldEnable_default = true;
//if(showFramesOn.contains(vi)){
auto polyscopeHandle_frameX = polyscopeHandle_deformedEdmeMesh
->addNodeVectorQuantity("FrameX", framesX);
polyscopeHandle_frameX->setVectorLengthScale(0.01);
polyscopeHandle_frameX->setVectorRadius(0.01);
polyscopeHandle_frameX->setVectorLengthScale(frameLength_default);
polyscopeHandle_frameX->setVectorRadius(frameRadius_default);
polyscopeHandle_frameX->setVectorColor(
/*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
auto polyscopeHandle_frameY = polyscopeHandle_deformedEdmeMesh
->addNodeVectorQuantity("FrameY", framesY);
polyscopeHandle_frameY->setVectorLengthScale(0.01);
polyscopeHandle_frameY->setVectorRadius(0.01);
polyscopeHandle_frameY->setVectorLengthScale(frameLength_default);
polyscopeHandle_frameY->setVectorRadius(frameRadius_default);
polyscopeHandle_frameY->setVectorColor(
/*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
auto polyscopeHandle_frameZ = polyscopeHandle_deformedEdmeMesh
->addNodeVectorQuantity("FrameZ", framesZ);
polyscopeHandle_frameZ->setVectorLengthScale(0.01);
polyscopeHandle_frameZ->setVectorRadius(0.01);
polyscopeHandle_frameZ->setVectorLengthScale(frameLength_default);
polyscopeHandle_frameZ->setVectorRadius(frameRadius_default);
polyscopeHandle_frameZ->setVectorColor(
/*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
//}
job->registerForDrawing(getLabel());
// auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel());
// if (!polyscopeHandle_initialMesh) {
// polyscopeHandle_initialMesh = mesh->registerForDrawing();
// }
// auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh
// ->addNodeVectorQuantity("FrameX", framesX_initial);
// polyscopeHandle_frameX_initial->setVectorLengthScale(frameLength_default);
// polyscopeHandle_frameX_initial->setVectorRadius(frameRadius_default);
// polyscopeHandle_frameX_initial->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
// auto polyscopeHandle_frameY_initial = polyscopeHandle_initialMesh
// ->addNodeVectorQuantity("FrameY", framesY_initial);
// polyscopeHandle_frameY_initial->setVectorLengthScale(frameLength_default);
// polyscopeHandle_frameY_initial->setVectorRadius(frameRadius_default);
// polyscopeHandle_frameY_initial->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
// auto polyscopeHandle_frameZ_initial = polyscopeHandle_initialMesh
// ->addNodeVectorQuantity("FrameZ", framesZ_initial);
// polyscopeHandle_frameZ_initial->setVectorLengthScale(frameLength_default);
// polyscopeHandle_frameZ_initial->setVectorRadius(frameRadius_default);
// polyscopeHandle_frameZ_initial->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
// //}
pJob->registerForDrawing(getLabel());
static bool wasExecuted = false;
if (!wasExecuted) {
std::function<void()> callback = [&]() {
@ -651,6 +732,7 @@ struct SimulationResults
PolyscopeInterface::addUserCallback(callback);
wasExecuted = true;
}
return polyscopeHandle_deformedEdmeMesh;
}
#endif
};

157
simulationmesh.cpp
View File

@ -1,5 +1,5 @@
#include "simulationmesh.hpp"
#include <wrap/nanoply/include/nanoplyWrapper.hpp>
//#include <wrap/nanoply/include/nanoplyWrapper.hpp>
SimulationMesh::SimulationMesh() {
elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(
@ -292,23 +292,23 @@ bool SimulationMesh::load(const string &plyFilename) {
this->Clear();
// assert(plyFileHasAllRequiredFields(plyFilename));
// Load the ply file
VCGEdgeMesh::PerEdgeAttributeHandle<CrossSectionType> handleBeamDimensions =
vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<
CrossSectionType>(*this, plyPropertyBeamDimensionsID);
VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial> handleBeamMaterial =
vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<ElementMaterial>(
*this, plyPropertyBeamMaterialID);
nanoply::NanoPlyWrapper<SimulationMesh>::CustomAttributeDescriptor
customAttrib;
customAttrib.GetMeshAttrib(plyFilename);
customAttrib.AddEdgeAttribDescriptor<CrossSectionType, double, 2>(
plyPropertyBeamDimensionsID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
/*FIXME: Since I allow CrossSectionType to take two types I should export the
* type as well such that that when loaded the correct type of cross section
* is used.
*/
customAttrib.AddEdgeAttribDescriptor<vcg::Point2d, double, 2>(
plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
// VCGEdgeMesh::PerEdgeAttributeHandle<CrossSectionType> handleBeamDimensions =
// vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<
// CrossSectionType>(*this, plyPropertyBeamDimensionsID);
// VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial> handleBeamMaterial =
// vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<ElementMaterial>(
// *this, plyPropertyBeamMaterialID);
// nanoply::NanoPlyWrapper<SimulationMesh>::CustomAttributeDescriptor
// customAttrib;
// customAttrib.GetMeshAttrib(plyFilename);
// customAttrib.AddEdgeAttribDescriptor<CrossSectionType, double, 2>(
// plyPropertyBeamDimensionsID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
// /*FIXME: Since I allow CrossSectionType to take two types I should export the
// * type as well such that that when loaded the correct type of cross section
// * is used.
// */
// customAttrib.AddEdgeAttribDescriptor<vcg::Point2d, double, 2>(
// plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
// VCGEdgeMesh::PerEdgeAttributeHandle<std::array<double, 6>>
// handleBeamProperties =
@ -322,21 +322,20 @@ bool SimulationMesh::load(const string &plyFilename) {
// customAttrib.AddEdgeAttribDescriptor<vcg::Point4f, float, 4>(
// plyPropertyBeamRigidityConstantsID, nanoply::NNP_LIST_INT8_FLOAT32,
// nullptr);
unsigned int mask = 0;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_MESHATTRIB;
if (nanoply::NanoPlyWrapper<SimulationMesh>::LoadModel(
plyFilename.c_str(), *this, mask, customAttrib) != 0) {
return false;
// unsigned int mask = 0;
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_MESHATTRIB;
if (!load(plyFilename)) {
return false;
}
elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(
*this, std::string("Elements"));
nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(
*this, std::string("Nodes"));
// elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(
// *this, std::string("Elements"));
// nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(
// *this, std::string("Nodes"));
vcg::tri::UpdateTopology<SimulationMesh>::VertexEdge(*this);
initializeNodes();
initializeElements();
@ -349,11 +348,11 @@ bool SimulationMesh::load(const string &plyFilename) {
// elements[ei].updateRigidity();
// }
// }
for (size_t ei = 0; ei < EN(); ei++) {
elements[ei].setDimensions(handleBeamDimensions[ei]);
elements[ei].setMaterial(handleBeamMaterial[ei]);
elements[ei].updateRigidity();
}
// for (size_t ei = 0; ei < EN(); ei++) {
// elements[ei].setDimensions(handleBeamDimensions[ei]);
// elements[ei].setMaterial(handleBeamMaterial[ei]);
// elements[ei].updateRigidity();
// }
bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001;
if (normalsAreAbsent) {
@ -375,28 +374,31 @@ bool SimulationMesh::save(const std::string &plyFilename)
if (filename.empty()) {
filename = std::filesystem::current_path().append(getLabel() + ".ply").string();
}
nanoply::NanoPlyWrapper<VCGEdgeMesh>::CustomAttributeDescriptor customAttrib;
customAttrib.GetMeshAttrib(filename);
// nanoply::NanoPlyWrapper<VCGEdgeMesh>::CustomAttributeDescriptor customAttrib;
// customAttrib.GetMeshAttrib(filename);
std::vector<CrossSectionType> dimensions = getBeamDimensions();
customAttrib.AddEdgeAttribDescriptor<CrossSectionType, double, 2>(plyPropertyBeamDimensionsID,
nanoply::NNP_LIST_INT8_FLOAT64,
dimensions.data());
std::vector<ElementMaterial> material = getBeamMaterial();
customAttrib.AddEdgeAttribDescriptor<vcg::Point2d, double, 2>(plyPropertyBeamMaterialID,
nanoply::NNP_LIST_INT8_FLOAT64,
material.data());
unsigned int mask = 0;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::SaveModel(filename.c_str(),
*this,
mask,
customAttrib,
false)
!= 1) {
// std::vector<CrossSectionType> dimensions = getBeamDimensions();
// customAttrib.AddEdgeAttribDescriptor<CrossSectionType, double, 2>(plyPropertyBeamDimensionsID,
// nanoply::NNP_LIST_INT8_FLOAT64,
// dimensions.data());
// std::vector<ElementMaterial> material = getBeamMaterial();
// customAttrib.AddEdgeAttribDescriptor<vcg::Point2d, double, 2>(plyPropertyBeamMaterialID,
// nanoply::NNP_LIST_INT8_FLOAT64,
// material.data());
// unsigned int mask = 0;
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
// if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::SaveModel(filename.c_str(),
// *this,
// mask,
// customAttrib,
// false)
// != 1) {
// return false;
// }
if (!save(plyFilename)) {
return false;
}
return true;
@ -449,31 +451,34 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1,
const VectorType cp = vector0 ^ vector1;
if (cp.dot(normalVector) < 0) {
return -angle;
return -angle;
}
return angle;
}
void Element::computeMaterialProperties(const ElementMaterial &material) {
this->G = material.youngsModulus / (2 * (1 + material.poissonsRatio));
G = material.youngsModulus / (2 * (1 + material.poissonsRatio));
}
void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A)
{
A = dimensions.b * dimensions.h;
}
void Element::computeDimensionsProperties(
const RectangularBeamDimensions &dimensions) {
assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions));
A = (dimensions.b * dimensions.h);
I2 = dimensions.b * std::pow(dimensions.h, 3) / 12;
I3 = dimensions.h * std::pow(dimensions.b, 3) / 12;
J = I2 + I3;
computeCrossSectionArea(dimensions, A);
computeMomentsOfInertia(dimensions, inertia);
}
void Element::computeDimensionsProperties(
const CylindricalBeamDimensions &dimensions) {
assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2));
I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64;
I3 = I2;
J = I2 + I3;
inertia.I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64;
inertia.I3 = inertia.I2;
inertia.J = inertia.I2 + inertia.I3;
}
void Element::setDimensions(const CrossSectionType &dimensions) {
@ -488,9 +493,23 @@ void Element::setMaterial(const ElementMaterial &material) {
updateRigidity();
}
double Element::getMass(const double &materialDensity)
{
const double beamVolume = A * length;
return beamVolume * materialDensity;
}
void Element::computeMomentsOfInertia(const RectangularBeamDimensions &dimensions,
Element::MomentsOfInertia &inertia)
{
inertia.I2 = dimensions.b * std::pow(dimensions.h, 3) / 12;
inertia.I3 = dimensions.h * std::pow(dimensions.b, 3) / 12;
inertia.J = inertia.I2 + inertia.I3;
}
void Element::updateRigidity() {
rigidity.axial = material.youngsModulus * A / initialLength;
rigidity.torsional = G * J / initialLength;
rigidity.firstBending = 2 * material.youngsModulus * I2 / initialLength;
rigidity.secondBending = 2 * material.youngsModulus * I3 / initialLength;
rigidity.torsional = G * inertia.J / initialLength;
rigidity.firstBending = 2 * material.youngsModulus * inertia.I2 / initialLength;
rigidity.secondBending = 2 * material.youngsModulus * inertia.I3 / initialLength;
}

16
simulationmesh.hpp
View File

@ -64,15 +64,23 @@ struct Element {
ElementMaterial material;
double G{0};
double A{0}; // cross sectional area
double I2{0}; // second moment of inertia
double I3{0}; // third moment of inertia
double J{0}; // torsional constant (polar moment of inertia)
void computeMaterialProperties(const ElementMaterial &material);
void computeDimensionsProperties(const RectangularBeamDimensions &dimensions);
void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions);
void setDimensions(const CrossSectionType &dimensions);
void setMaterial(const ElementMaterial &material);
double getMass(const double &matrialDensity);
struct MomentsOfInertia
{
double I2{0}; // second moment of inertia
double I3{0}; // third moment of inertia
double J{0}; // torsional constant (polar moment of inertia)
} inertia;
static void computeMomentsOfInertia(const RectangularBeamDimensions &dimensions,
MomentsOfInertia &inertia);
struct LocalFrame {
VectorType t1;
@ -127,6 +135,8 @@ struct Element {
};
RotationalDisplacements rotationalDisplacements_j;
RotationalDisplacements rotationalDisplacements_jplus1;
static void computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A);
};
struct Node {

13
simulationmodel.hpp Normal file
View File

@ -0,0 +1,13 @@
#ifndef SIMULATIONMODEL_HPP
#define SIMULATIONMODEL_HPP
#include "simulation_structs.hpp"
class SimulationModel
{
public:
virtual SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &simulationJob)
= 0;
};
#endif // SIMULATIONMODEL_HPP

77
topologyenumerator.cpp
View File

@ -152,6 +152,9 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
}
}
const std::unordered_set<VertexIndex> interfaceNodes = patternGeometryAllEdges.getInterfaceNodes(
numberOfNodesPerSlot);
// assert(validEdges.size() == allPossibleEdges.size() -
// coincideEdges.size() -
// duplicateEdges.size());
@ -191,7 +194,8 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
perEdgeResultPath,
patternGeometryAllEdges.getVertices(),
intersectingEdges,
validEdges);
validEdges,
interfaceNodes);
statistics.print(setupString, perEdgeResultPath);
statistics.reset();
}
@ -211,7 +215,8 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
perEdgeResultPath,
patternGeometryAllEdges.getVertices(),
intersectingEdges,
validEdges);
validEdges,
interfaceNodes);
statistics.print(setupString, perEdgeResultPath);
}
}
@ -395,7 +400,8 @@ void TopologyEnumerator::computeValidPatterns(
const std::filesystem::path &resultsPath,
const std::vector<vcg::Point3d> &allVertices,
const std::unordered_map<size_t, std::unordered_set<size_t>> &intersectingEdges,
const std::vector<vcg::Point2i> &validEdges)
const std::vector<vcg::Point2i> &validEdges,
const std::unordered_set<VertexIndex> &interfaceNodes)
{
assert(numberOfNodesPerSlot.size() == 7);
// Iterate over all patterns which have numberOfDesiredEdges edges from
@ -435,7 +441,7 @@ void TopologyEnumerator::computeValidPatterns(
// std::string previousPatternBinaryRepresentation(validEdges.size(),'0');
size_t patternIndex = 0;
bool validPatternsExist = false;
const bool exportTilledPattern = debugIsOn;
const bool exportTilledPattern = true;
const bool saveCompressedFormat = false;
do {
patternIndex++;
@ -459,6 +465,23 @@ void TopologyEnumerator::computeValidPatterns(
patternGeometry.add(allVertices, patternEdges);
patternGeometry.setLabel(patternName);
// Check if pattern contains intersecting edges
const bool isInterfaceConnected = patternGeometry.isInterfaceConnected(interfaceNodes);
// Export the tiled ply file if it contains intersecting edges
if (!isInterfaceConnected) {
// create the tiled geometry of the pattern
statistics.numberOfPatternViolatingInterfaceEnforcement++;
if (debugIsOn) {
if (savePlyFiles) {
exportPattern(std::filesystem::path(resultsPath).append("InterfaceEnforcement"),
patternGeometry,
exportTilledPattern);
}
} else {
continue; // should be uncommented in order to improve performance
}
}
// Check if pattern contains intersecting edges
const bool patternContainsIntersectingEdges
= patternGeometry.hasIntersectingEdges(patternBinaryRepresentation, intersectingEdges);
@ -477,6 +500,7 @@ void TopologyEnumerator::computeValidPatterns(
}
}
// const bool shouldBreak = numberOfDesiredEdges == 4 && patternIndex == 53;
const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold
= patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 15);
if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) {
@ -512,7 +536,7 @@ void TopologyEnumerator::computeValidPatterns(
const bool tiledPatternHasDanglingEdges = patternGeometry.hasDanglingEdges(
numberOfNodesPerSlot); // marks the nodes with valence>=1
// Create the tiled geometry of the pattern
const bool hasFloatingComponents = !patternGeometry.isFullyConnectedWhenTiled();
const bool hasFloatingComponents = !patternGeometry.isFullyConnectedWhenFanned();
PatternGeometry fanPatternGeometry = PatternGeometry::createFan(patternGeometry);
const int interfaceNodeVi = 3;
@ -559,10 +583,19 @@ void TopologyEnumerator::computeValidPatterns(
std::vector<vertex_t> articulationPoints;
boost::articulation_points(fanPatternGraph, std::back_inserter(articulationPoints));
const bool hasArticulationPoints = !articulationPoints.empty();
// if (hasArticulationPoints /*&& !patternContainsIntersectingEdges
// && !tiledPatternHasDanglingEdges && !hasFloatingComponents
// && !tiledPatternHasNodeWithValenceGreaterThanDesired
// && !tiledPatternHasEdgesWithAngleSmallerThanThreshold*/) {
// if (!hasArticulationPoints && tiledPatternHasDanglingEdges) {
// PatternGeometry tilledPatternGeometry = PatternGeometry::createTile(patternGeometry);
// tilledPatternGeometry.updateEigenEdgeAndVertices();
// tilledPatternGeometry.registerForDrawing();
// // ->addNodeColorQuantity("de_noAp_tilled", fanVertexColors)
// // ->setEnabled(true);
// polyscope::show();
// tilledPatternGeometry.unregister();
// }
// if (hasArticulationPoints && !tiledPatternHasDanglingEdges/*&& !patternContainsIntersectingEdges
// && !hasFloatingComponents
// && !tiledPatternHasNodeWithValenceGreaterThanDesired
// && !tiledPatternHasEdgesWithAngleSmallerThanThreshold*/) {
// for (PatternGeometry::VertexType &v : patternGeometry.vert) {
// v.C() = vcg::Color4b::Yellow;
// }
@ -582,12 +615,12 @@ void TopologyEnumerator::computeValidPatterns(
// const auto vi = tilledPatternGeometry.getIndex(v);
// fanVertexColors[vi] = vColor;
// }
// // tilledPatternGeometry.updateEigenEdgeAndVertices();
// // tilledPatternGeometry.registerForDrawing()
// // ->addNodeColorQuantity("ap_tilled", fanVertexColors)
// // ->setEnabled(true);
// // polyscope::show();
// // tilledPatternGeometry.unregister();
// tilledPatternGeometry.updateEigenEdgeAndVertices();
// tilledPatternGeometry.registerForDrawing()
// ->addNodeColorQuantity("ap_tilled", fanVertexColors)
// ->setEnabled(true);
// polyscope::show();
// tilledPatternGeometry.unregister();
// }
// duplicated here
// Check dangling edges with vcg method
@ -687,10 +720,20 @@ void TopologyEnumerator::computeValidPatterns(
}
const bool isValidPattern = !patternContainsIntersectingEdges
&& !tiledPatternHasDanglingEdges && !hasFloatingComponents
&& !hasArticulationPoints
&& isInterfaceConnected
/*&& !tiledPatternHasDanglingEdges*/
&& !hasFloatingComponents && !hasArticulationPoints
&& !tiledPatternHasNodeWithValenceGreaterThanDesired
&& !tiledPatternHasEdgesWithAngleSmallerThanThreshold;
// if (!hasArticulationPoints && !patternContainsIntersectingEdges
// && !tiledPatternHasDanglingEdges && !hasFloatingComponents
// && !tiledPatternHasNodeWithValenceGreaterThanDesired
// && tiledPatternHasEdgesWithAngleSmallerThanThreshold && numberOfDesiredEdges > 4) {
// std::cout << "Pattern found:" << patternName << std::endl;
// patternGeometry.registerForDrawing();
// polyscope::show();
// patternGeometry.unregister();
// }
if (isValidPattern) {
// if(patternName=='2055'){
// PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(

26
topologyenumerator.hpp
View File

@ -64,25 +64,27 @@ class TopologyEnumerator {
size_t numberOfPatternsWithArticulationPoints{0};
size_t numberOfPatternsViolatingAngleThreshold{0};
size_t numberOfPatternsViolatingValenceThreshold{0};
size_t numberOfPatternViolatingInterfaceEnforcement{0};
size_t numberOfValidPatterns{0};
nlohmann::json convertToJson() const
{
nlohmann::json json;
json["numPossibleEdges"] = numberOfPossibleEdges;
json["numCoincideEdges"] = numberOfCoincideEdges;
json["numDuplicateEdges"] = numberOfDuplicateEdges;
json["numValidEdges"] = numberOfValidEdges;
json["numIntersectingEdgePairs"] = numberOfIntersectingEdgePairs;
// json["numPossibleEdges"] = numberOfPossibleEdges;
// json["numCoincideEdges"] = numberOfCoincideEdges;
// json["numDuplicateEdges"] = numberOfDuplicateEdges;
// json["numValidEdges"] = numberOfValidEdges;
// json["numIntersectingEdgePairs"] = numberOfIntersectingEdgePairs;
json["numPatterns"] = numberOfPatterns;
// json["numIntersectingEdgesOverAllPatterns"] =
// numberOfIntersectingEdgesOverAllPatterns;
json["numPatternsWithIntersectingEdges"] = numberOfPatternsWithIntersectingEdges;
json["numViolatingInterfaceEnforcement"] = numberOfPatternViolatingInterfaceEnforcement;
json["numPatternsWithNotASingleCC"] = numberOfPatternsWithMoreThanASingleCC;
json["numPatternsWithDangling"] = numberOfPatternsWithADanglingEdgeOrNode;
json["numPatternsWithArticulationPoints"] = numberOfPatternsWithArticulationPoints;
json["numValidPatterns"] = numberOfValidPatterns;
json["numPatternsWithDangling"] = numberOfPatternsWithADanglingEdgeOrNode;
json["violatingAngle"] = numberOfPatternsViolatingAngleThreshold;
json["violatingValence"] = numberOfPatternsViolatingValenceThreshold;
json["numValidPatterns"] = numberOfValidPatterns;
return json;
}
@ -107,9 +109,10 @@ class TopologyEnumerator {
// std::endl;
std::cout << numberOfPatternsWithIntersectingEdges
<< " patterns found with intersecting edges" << std::endl;
std::cout << numberOfPatternViolatingInterfaceEnforcement
<< " patterns found which violate interface enforcement" << std::endl;
std::cout << numberOfPatternsWithMoreThanASingleCC
<< " patterns found with more than one connected components"
<< std::endl;
<< " patterns found with more than one connected components" << std::endl;
std::cout << numberOfPatternsWithADanglingEdgeOrNode
<< " patterns found with a dangling node or edge" << std::endl;
std::cout << numberOfPatternsWithArticulationPoints
@ -141,6 +144,8 @@ class TopologyEnumerator {
numberOfIntersectingEdgePairs = 0;
numberOfPatterns = 0;
numberOfPatternsWithIntersectingEdges = 0;
numberOfPatternViolatingInterfaceEnforcement = 0;
;
numberOfPatternsWithMoreThanASingleCC = 0;
numberOfPatternsWithADanglingEdgeOrNode = 0;
numberOfPatternsWithArticulationPoints = 0;
@ -195,7 +200,8 @@ private:
const std::filesystem::path &resultsPath,
const std::vector<vcg::Point3d> &vertices,
const std::unordered_map<size_t, std::unordered_set<size_t>> &intersectingEdges,
const std::vector<vcg::Point2i> &validEdges);
const std::vector<vcg::Point2i> &validEdges,
const std::unordered_set<VertexIndex> &interfaceNodes);
void exportPattern(const std::filesystem::path &saveToPath,
PatternGeometry &patternGeometry,
const bool saveTilledPattern) const;

248
trianglepatterngeometry.cpp
View File

@ -178,10 +178,11 @@ void PatternGeometry::add(const std::vector<vcg::Point2i> &edges) {
}
void PatternGeometry::add(const std::vector<vcg::Point3d> &vertices,
const std::vector<vcg::Point2i> &edges) {
add(vertices);
add(edges);
updateEigenEdgeAndVertices();
const std::vector<vcg::Point2i> &edges)
{
add(vertices);
add(edges);
updateEigenEdgeAndVertices();
}
void PatternGeometry::add(const std::vector<size_t> &numberOfNodesPerSlot,
@ -375,9 +376,9 @@ bool PatternGeometry::hasAngleSmallerThanThreshold(const std::vector<size_t> &nu
// const CoordType &correspondingVertexPosition = vert[correspondingNode[vi]].cP();
vcg::Matrix33d R;
if (slotIndex == 3) {
R = vcg::RotationMatrix(vcg::Point3d(0, 0, 1), vcg::math::ToRad(-60.0));
R = vcg::RotationMatrix(vcg::Point3d(0, 0, 1), vcg::math::ToRad(-360.0 / fanSize));
} else {
R = vcg::RotationMatrix(vcg::Point3d(0, 0, 1), vcg::math::ToRad(60.0));
R = vcg::RotationMatrix(vcg::Point3d(0, 0, 1), vcg::math::ToRad(360.0 / fanSize));
}
// const CoordType &correspondingVertexPosition_rotated = vert[vi].cP();
std::transform(
@ -452,6 +453,7 @@ bool PatternGeometry::hasAngleSmallerThanThreshold(const std::vector<size_t> &nu
}
return false;
}
bool PatternGeometry::hasValenceGreaterThan(const std::vector<size_t> &numberOfNodesPerSlot,
const size_t &valenceThreshold)
{
@ -561,111 +563,111 @@ void PatternGeometry::constructCorrespondingNodeMap(const std::vector<size_t> &n
}
}
bool PatternGeometry::isFullyConnectedWhenTiled() {
assert(vn != 0 /* && !correspondingNode.empty()*/);
// TrianglePatternGeometry copyOfPattern(*this);
bool PatternGeometry::isFullyConnectedWhenFanned()
{
assert(vn != 0 /* && !correspondingNode.empty()*/);
// TrianglePatternGeometry copyOfPattern(*this);
// If bottom interface nodes have a valence of zero there definetely more than
// a single cc
bool bottomInterfaceIsConnected = false;
for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
if (nodeToSlotMap[nodeIndex] == 1 || nodeToSlotMap[nodeIndex] == 4
|| nodeToSlotMap[nodeIndex] == 2) {
std::vector<PatternGeometry::EdgeType *> connectedEdges;
vcg::edge::VEStarVE(&vert[nodeIndex], connectedEdges);
const size_t nodeValence = connectedEdges.size();
if (nodeValence != 0) {
bottomInterfaceIsConnected = true;
break;
}
}
}
if (!bottomInterfaceIsConnected) {
return false;
}
// // If bottom interface nodes have a valence of zero there definetely more than
// // a single cc
// bool bottomInterfaceIsConnected = false;
// for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
// if (nodeToSlotMap[nodeIndex] == 1 || nodeToSlotMap[nodeIndex] == 4
// || nodeToSlotMap[nodeIndex] == 2) {
// std::vector<PatternGeometry::EdgeType *> connectedEdges;
// vcg::edge::VEStarVE(&vert[nodeIndex], connectedEdges);
// const size_t nodeValence = connectedEdges.size();
// if (nodeValence != 0) {
// bottomInterfaceIsConnected = true;
// break;
// }
// }
// }
// if (!bottomInterfaceIsConnected) {
// return false;
// }
PatternGeometry fanedPattern = createFan(*this);
vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(fanedPattern, 0.000000005);
vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(fanedPattern);
vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(fanedPattern);
vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(fanedPattern);
std::vector<std::pair<int, PatternGeometry::EdgePointer>> eCC;
vcg::tri::Clean<PatternGeometry>::edgeMeshConnectedComponents(fanedPattern,
eCC);
PatternGeometry fanedPattern = createFan(*this);
vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(fanedPattern, 0.000000005);
vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(fanedPattern);
vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(fanedPattern);
vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(fanedPattern);
std::vector<std::pair<int, PatternGeometry::EdgePointer>> eCC;
vcg::tri::Clean<PatternGeometry>::edgeMeshConnectedComponents(fanedPattern, eCC);
const bool sideInterfaceIsConnected = 1 == eCC.size();
if (!sideInterfaceIsConnected) {
return false;
}
const bool sideInterfaceIsConnected = 1 == eCC.size();
if (!sideInterfaceIsConnected) {
return false;
}
// for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
// if (nodeTiledValence[nodeIndex] == 0)
// continue;
// const size_t slotIndex = nodeSlot[nodeIndex];
// // connect nodes belonging to first or third slot(nodes on the first
// // triangle edge)
// // if (nodeTiledValence[nodeIndex] == 0 && slotIndex == 3) {
// // continue;
// // }
// if (slotIndex == 1 || slotIndex == 3) {
// assert(correspondingNode.count(nodeIndex) != 0);
// const size_t correspondingNodeIndex =
// correspondingNode[nodeIndex];
// std::vector<TrianglePatternGeometry::EdgeType *> starEdges;
// vcg::edge::VEStarVE(&vert[nodeIndex], starEdges);
// bool containsEdge = false;
// for (TrianglePatternGeometry::EdgeType *e : starEdges) {
// assert(vcg::tri::Index(*this, e->V(0)) == nodeIndex ||
// vcg::tri::Index(*this, e->V(1)) == nodeIndex);
// if (vcg::tri::Index(*this, e->V(0)) == nodeIndex) {
// if (vcg::tri::Index(*this, e->V(1)) == correspondingNodeIndex) {
// containsEdge = true;
// break;
// }
// } else if (vcg::tri::Index(*this, e->V(1)) == nodeIndex) {
// if (vcg::tri::Index(*this, e->V(0)) == correspondingNodeIndex) {
// containsEdge = true;
// break;
// }
// }
// }
// if (!containsEdge) {
// vcg::tri::Allocator<TrianglePatternGeometry>::AddEdge(
// copyOfPattern, nodeIndex, correspondingNodeIndex);
// }
// } else if (slotIndex == 2 || slotIndex == 5) {
// assert(correspondingNode.count(nodeIndex) != 0);
// } else {
// assert(correspondingNode.count(nodeIndex) == 0);
// }
// }
// for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
// if (nodeTiledValence[nodeIndex] == 0)
// continue;
// const size_t slotIndex = nodeSlot[nodeIndex];
// // connect nodes belonging to first or third slot(nodes on the first
// // triangle edge)
// // if (nodeTiledValence[nodeIndex] == 0 && slotIndex == 3) {
// // continue;
// // }
// if (slotIndex == 1 || slotIndex == 3) {
// assert(correspondingNode.count(nodeIndex) != 0);
// const size_t correspondingNodeIndex =
// correspondingNode[nodeIndex];
// std::vector<TrianglePatternGeometry::EdgeType *> starEdges;
// vcg::edge::VEStarVE(&vert[nodeIndex], starEdges);
// bool containsEdge = false;
// for (TrianglePatternGeometry::EdgeType *e : starEdges) {
// assert(vcg::tri::Index(*this, e->V(0)) == nodeIndex ||
// vcg::tri::Index(*this, e->V(1)) == nodeIndex);
// if (vcg::tri::Index(*this, e->V(0)) == nodeIndex) {
// if (vcg::tri::Index(*this, e->V(1)) == correspondingNodeIndex) {
// containsEdge = true;
// break;
// }
// } else if (vcg::tri::Index(*this, e->V(1)) == nodeIndex) {
// if (vcg::tri::Index(*this, e->V(0)) == correspondingNodeIndex) {
// containsEdge = true;
// break;
// }
// }
// }
// if (!containsEdge) {
// vcg::tri::Allocator<TrianglePatternGeometry>::AddEdge(
// copyOfPattern, nodeIndex, correspondingNodeIndex);
// }
// } else if (slotIndex == 2 || slotIndex == 5) {
// assert(correspondingNode.count(nodeIndex) != 0);
// } else {
// assert(correspondingNode.count(nodeIndex) == 0);
// }
// }
// std::vector<std::pair<int, TrianglePatternGeometry::EdgePointer>> eCC;
// vcg::tri::Clean<TrianglePatternGeometry>::edgeMeshConnectedComponents(
// copyOfPattern, eCC);
// size_t numberOfCC_edgeBased = eCC.size();
// size_t numberOfCC_vertexBased = numberOfCC_edgeBased;
// if (numberOfCC_edgeBased == 1) {
// vcg::tri::UpdateTopology<TrianglePatternGeometry>::VertexEdge(
// copyOfPattern);
// vcg::tri::UpdateTopology<TrianglePatternGeometry>::EdgeEdge(copyOfPattern);
// vcg::tri::UpdateFlags<TrianglePatternGeometry>::VertexSetV(copyOfPattern);
// vcg::tri::UpdateFlags<TrianglePatternGeometry>::VertexClear(copyOfPattern);
// for (size_t ei = 0; ei < copyOfPattern.EN(); ei++) {
// copyOfPattern.edge[ei].V(0)->SetV();
// copyOfPattern.edge[ei].V(1)->SetV();
// }
// std::vector<std::pair<int, TrianglePatternGeometry::EdgePointer>> eCC;
// vcg::tri::Clean<TrianglePatternGeometry>::edgeMeshConnectedComponents(
// copyOfPattern, eCC);
// size_t numberOfCC_edgeBased = eCC.size();
// size_t numberOfCC_vertexBased = numberOfCC_edgeBased;
// if (numberOfCC_edgeBased == 1) {
// vcg::tri::UpdateTopology<TrianglePatternGeometry>::VertexEdge(
// copyOfPattern);
// vcg::tri::UpdateTopology<TrianglePatternGeometry>::EdgeEdge(copyOfPattern);
// vcg::tri::UpdateFlags<TrianglePatternGeometry>::VertexSetV(copyOfPattern);
// vcg::tri::UpdateFlags<TrianglePatternGeometry>::VertexClear(copyOfPattern);
// for (size_t ei = 0; ei < copyOfPattern.EN(); ei++) {
// copyOfPattern.edge[ei].V(0)->SetV();
// copyOfPattern.edge[ei].V(1)->SetV();
// }
// for (size_t vi = 0; vi < copyOfPattern.VN(); vi++) {
// if (!copyOfPattern.vert[vi].IsV()) {
// numberOfCC_vertexBased++;
// }
// }
// return numberOfCC_vertexBased;
// }
// for (size_t vi = 0; vi < copyOfPattern.VN(); vi++) {
// if (!copyOfPattern.vert[vi].IsV()) {
// numberOfCC_vertexBased++;
// }
// }
// return numberOfCC_vertexBased;
// }
// return numberOfCC_edgeBased == 1; // TODO: not good
return true;
// return numberOfCC_edgeBased == 1; // TODO: not good
return true;
}
bool PatternGeometry::hasIntersectingEdges(
@ -1062,7 +1064,7 @@ void PatternGeometry::createFan(const std::vector<int> &connectToNeighborsVi, co
for (int rotationCounter = 1; rotationCounter < fanSize; rotationCounter++) {
vcg::Matrix44d R;
auto rotationAxis = vcg::Point3d(0, 0, 1);
R.SetRotateDeg(360 / fanSize, rotationAxis);
R.SetRotateDeg(360.0 / fanSize, rotationAxis);
vcg::tri::UpdatePosition<PatternGeometry>::Matrix(rotatedPattern, R);
vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(*this, rotatedPattern);
//Add edges for connecting the desired vertices
@ -1097,3 +1099,39 @@ double PatternGeometry::getBaseTriangleHeight() const
{
return baseTriangleHeight;
}
bool PatternGeometry::isInterfaceConnected(const std::unordered_set<VertexIndex> &interfaceNodes)
{
std::unordered_set<VertexIndex> interfaceNodesConnected;
for (int ei = 0; ei < EN(); ei++) {
const int evi0 = getIndex(edge[ei].cV(0));
const int evi1 = getIndex(edge[ei].cV(1));
if (interfaceNodes.contains(evi0) && !interfaceNodes.contains(evi1)) {
interfaceNodesConnected.insert(evi0);
} else if (!interfaceNodes.contains(evi0) && interfaceNodes.contains(evi1)) {
interfaceNodesConnected.insert(evi1);
}
}
if (interfaceNodesConnected.size() != interfaceNodes.size()) {
return false;
}
return true;
}
std::unordered_set<VertexIndex> PatternGeometry::getInterfaceNodes(
const std::vector<size_t> &numberOfNodesPerSlot)
{
if (nodeToSlotMap.empty()) {
FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlotMap);
}
std::unordered_set<VertexIndex> interfaceNodes;
for (std::pair<VertexIndex, size_t> viSlotPair : nodeToSlotMap) {
if (viSlotPair.second == 4) {
interfaceNodes.insert(viSlotPair.first);
}
}
return interfaceNodes;
}

7
trianglepatterngeometry.hpp
View File

@ -49,8 +49,7 @@ private:
bool load(const std::string &plyFilename) override;
void add(const std::vector<vcg::Point3d> &vertices);
void add(const std::vector<vcg::Point2i> &edges);
void add(const std::vector<vcg::Point3d> &vertices,
const std::vector<vcg::Point2i> &edges);
void add(const std::vector<vcg::Point3d> &vertices, const std::vector<vcg::Point2i> &edges);
void add(const std::vector<size_t> &numberOfNodesPerSlot,
const std::vector<vcg::Point2i> &edges);
static std::vector<vcg::Point3d>
@ -71,7 +70,7 @@ private:
assert(n >= m);
return tgamma(n + 1) / (tgamma(m + 1) * tgamma(n - m + 1));
}
bool isFullyConnectedWhenTiled();
bool isFullyConnectedWhenFanned();
bool hasIntersectingEdges(
const std::string &patternBinaryRepresentation,
@ -123,6 +122,8 @@ private:
const std::vector<int> &perSurfaceFacePatternIndices,
std::vector<size_t> &tileIntoEdgesToTiledVi,
std::vector<std::vector<size_t>> &perPatternIndexTiledPatternEdgeIndex);
std::unordered_set<VertexIndex> getInterfaceNodes(const std::vector<size_t> &numberOfNodesPerSlot);
bool isInterfaceConnected(const std::unordered_set<VertexIndex> &interfaceNodes);
};
#endif // FLATPATTERNGEOMETRY_HPP

14
utilities.hpp
View File

@ -257,6 +257,7 @@ inline void init()
polyscope::view::upDir = polyscope::view::UpDir::ZUp;
polyscope::state::userCallback = &mainCallback;
// polyscope::options::autocenterStructures = true;
}
using PolyscopeLabel = std::string;
inline std::pair<PolyscopeLabel, size_t> getSelection()
@ -346,4 +347,17 @@ inline size_t computeHashOrdered(const std::vector<int> &v)
return std::hash<std::string>{}(elementsString);
}
//inline std::filesystem::path getFilepathWithExtension(const std::filesystem::path &folderPath,
// const std::string &extension)
//{
// for (const std::filesystem::directory_entry &dirEntry :
// std::filesystem::directory_iterator(folderPath)) {
// if (dirEntry.is_regular_file() && std::filesystem::path(dirEntry).extension() == extension) {
// return std::filesystem::path(dirEntry);
// }
// }
// return "";
//}
#endif // UTILITIES_H

37
vcgtrimesh.cpp
View File

@ -3,17 +3,18 @@
#include "wrap/io_trimesh/import_obj.h"
#include "wrap/io_trimesh/import_off.h"
#include <filesystem>
#include <wrap/io_trimesh/export.h>
#include <wrap/io_trimesh/import.h>
#include <wrap/nanoply/include/nanoplyWrapper.hpp>
//#include <wrap/nanoply/include/nanoplyWrapper.hpp>
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;
mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
mask |= vcg::tri::io::Mask::IOM_VERTCOLOR;
mask |= vcg::tri::io::Mask::IOM_EDGEINDEX;
mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
// if (nanoply::NanoPlyWrapper<VCGTriMesh>::LoadModel(
// std::filesystem::absolute(filename).string().c_str(), *this, mask)
@ -79,12 +80,11 @@ bool VCGTriMesh::save(const std::string plyFilename)
assert(std::filesystem::path(plyFilename).extension() == ".ply");
// Load the ply file
unsigned int mask = 0;
mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTCOORD;
mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_VERTCOLOR;
mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_FACEINDEX;
mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_FACENORMAL;
if (nanoply::NanoPlyWrapper<VCGTriMesh>::SaveModel(plyFilename.c_str(), *this, mask, false)
!= 0) {
mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
mask |= vcg::tri::io::Mask::IOM_VERTCOLOR;
mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
mask |= vcg::tri::io::Mask::IOM_FACENORMAL;
if (vcg::tri::io::Exporter<VCGTriMesh>::Save(*this, plyFilename.c_str(), mask) != 0) {
return false;
}
return true;
@ -117,7 +117,7 @@ polyscope::SurfaceMesh *VCGTriMesh::registerForDrawing(
const glm::vec3 desiredColor_glm(desiredColor.value()[0],
desiredColor.value()[1],
desiredColor.value()[2]);
polyscopeHandle_mesh->setSurfaceColor(desiredColor_glm);
polyscopeHandle_mesh->setSurfaceColor(glm::normalize(desiredColor_glm));
}
return polyscopeHandle_mesh;
}
@ -146,3 +146,14 @@ VCGTriMesh::VCGTriMesh(const std::string &filename)
vcg::tri::UpdateTopology<VCGTriMesh>::VertexFace(*this);
vcg::tri::UpdateNormal<VCGTriMesh>::PerVertexNormalized(*this);
}
void VCGTriMesh::moveToCenter()
{
CoordType centerOfMass(0, 0, 0);
for (int vi = 0; vi < VN(); vi++) {
centerOfMass += vert[vi].cP();
}
centerOfMass /= VN();
vcg::tri::UpdatePosition<VCGTriMesh>::Translate(*this, -centerOfMass);
}

1
vcgtrimesh.hpp
View File

@ -54,6 +54,7 @@ public:
#endif
Eigen::MatrixX2i getEdges() const;
void updateEigenEdgeAndVertices();
void moveToCenter();
};
#endif // VCGTRIMESH_HPP