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

4 parameters optimization

This commit is contained in:
Iason 2020-12-16 21:31:58 +02:00
parent b0ba64d499
commit b32e7ff8a9
3 changed files with 256 additions and 205 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

156
src/main.cpp
View File

@ -43,6 +43,7 @@ int main(int argc, char *argv[]) {
// pattern.savePly("fannedValid.ply");
// Create reduced models
// FormFinder::runUnitTests();
const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
std::vector<vcg::Point2i> singleBarReducedModelEdges{vcg::Point2i(0, 3)};
FlatPattern singleBarReducedModel(numberOfNodesPerSlot,
@ -63,6 +64,14 @@ int main(int argc, char *argv[]) {
&CWReducedModel, &CCWReducedModel};
ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
// FlatPattern pattern("/home/iason/Models/TestSet_validPatterns/59.ply");
// optimizer.initialize(pattern, *reducedModels[1], {});
// optimizer.setInitialGuess({0.10073746137216238, 1.9611251949927333,
// 7.5940127902102867, 8.0260424171674654});
// Eigen::VectorXd optimalParameters = optimizer.optimize();
// optimizer.setInitialGuess({2, 2, 2, 2, 2});
std::string fullPatternsTestSetDirectory =
"/home/iason/Models/TestSet_validPatterns";
for (const auto &entry :
@ -75,13 +84,14 @@ int main(int argc, char *argv[]) {
// std::filesystem::path(
// "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/"
// "1v_0v_2e_1e_1c_6fan/2/Valid/33.ply");
std::filesystem::path(
"/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/"
"1v_0v_2e_1e_1c_6fan/3/Valid/431.ply");
// std::filesystem::path(
// "/home/iason/Models/TestSet_validPatterns/865.ply");
// "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/"
// "1v_0v_2e_1e_1c_6fan/3/Valid/431.ply");
// std::filesystem::path(
// "/home/iason/Models/TestSet_validPatterns/59.ply");
entry.path();
const auto filepathString = filepath.string();
std::cout << "Full pattern:" << filepathString << std::endl;
// Use only the base triangle version
const std::string tiledSuffix = "_tiled.ply";
if (filepathString.compare(filepathString.size() - tiledSuffix.size(),
@ -91,18 +101,23 @@ int main(int argc, char *argv[]) {
FlatPattern pattern(filepathString);
pattern.setLabel(filepath.stem().string());
std::cout << "Testing Pattern:" << filepathString << std::endl;
for (FlatPattern *pReducedModel : reducedModels) {
// pReducedModel = reducedModels[1];
// for (int reducedPatternIndex = 0;
// reducedPatternIndex < reducedModels.size();
// reducedPatternIndex++) {
// FlatPattern *pReducedModel = reducedModels[reducedPatternIndex];
std::unordered_set<size_t> optimizationExcludedEi;
if (pReducedModel !=
reducedModels[0]) { // assumes that the singleBar reduced model is
// the
// first in the reducedModels vector
optimizationExcludedEi.insert(0);
}
optimizer.initialize(pattern, *pReducedModel, optimizationExcludedEi);
// if (pReducedModel !=
// reducedModels[0]) { // assumes that the singleBar reduced model
// // is
// // the
// // first in the reducedModels vector
// optimizationExcludedEi.insert(0);
// }
std::cout << "Full pattern:" << filepathString << std::endl;
std::cout << "Reduced pattern:" << 1 << std::endl;
optimizer.initialize(pattern, *reducedModels[1], optimizationExcludedEi);
Eigen::VectorXd optimalParameters = optimizer.optimize();
}
// }
}
// // Full model simulation
@ -170,62 +185,65 @@ int main(int argc, char *argv[]) {
// stiffnessFactor *= 1.5;
// }
// Beam
// // Beam
// // registerWorldAxes();
// VCGEdgeMesh mesh;
// // mesh.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply");
// mesh.loadFromPly("/home/iason/Coding/Projects/Approximating shapes with
// flat "
// "patterns/RodModelOptimizationForPatterns/build/Debug/"
// "Fanned_Single bar reduced model.ply");
// // vcg::Matrix44d R;
// // R.SetRotateDeg(45, {0, 0, 1});
// // vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(mesh, R);
// // mesh.updateEigenEdgeAndVertices();
// FormFinder formFinder;
// formFinder.runUnitTests();
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
// fixedVertices; fixedVertices[1] = std::unordered_set<DoFType>{2};
// fixedVertices[2] = std::unordered_set<DoFType>{2};
// fixedVertices[3] = std::unordered_set<DoFType>{1, 2};
// fixedVertices[4] = std::unordered_set<DoFType>{2};
// fixedVertices[5] = std::unordered_set<DoFType>{2};
// fixedVertices[6] = std::unordered_set<DoFType>{0, 1, 2};
// // Forced displacements
// std::unordered_map<size_t, Eigen::Vector3d> nodalForcedDisplacements;
// // nodalForcedDisplacements[10] = Eigen::Vector3d(0, 0.1, 0.1);
// std::unordered_map<size_t, VectorType> nodalForcedNormal;
// CoordType v14 = (mesh.vert[4].cP() - mesh.vert[0].cP()).Normalize();
// CoordType v25 = (mesh.vert[5].cP() - mesh.vert[0].cP()).Normalize();
// CoordType v36 = (mesh.vert[6].cP() - mesh.vert[0].cP()).Normalize();
// const double forceMagnitude = 0.4;
// std::unordered_map<VertexIndex, Vector6d> nodalForces{
// // {4, Vector6d({0, 0, 0, v14[0], v14[1], 0}) * forceMagnitude},
// // {1, Vector6d({0, 0, 0, -v14[0], -v14[1], 0}) *
// forceMagnitude},
// // {5, Vector6d({0, 0, 0, v25[0], v25[1], 0}) * forceMagnitude},
// // {2, Vector6d({0, 0, 0, -v25[0], -v25[1], 0}) *
// forceMagnitude},
// // {6, Vector6d({0, 0, 0, v36[0], v36[1], 0}) * forceMagnitude},
// {3, Vector6d({0, 0, 0, -v36[0], -v36[1], 0}) * forceMagnitude}};
// // nodalForcedNormal[10] = v;
// // nodalForcedNormal[0] = -v;
// // fixedVertices[10] = {0, 1};
// SimulationJob beamSimulationJob{std::make_shared<SimulationMesh>(mesh),
// fixedVertices,
// nodalForces,
// {},
// {}};
// beamSimulationJob.mesh->setBeamMaterial(0.3, 1);
// beamSimulationJob.mesh->setBeamCrossSection(
// RectangularBeamDimensions{0.002, 0.002});
// beamSimulationJob.registerForDrawing();
// registerWorldAxes();
VCGEdgeMesh mesh;
// mesh.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply");
mesh.loadFromPly("/home/iason/Coding/Projects/Approximating shapes with flat "
"patterns/RodModelOptimizationForPatterns/build/Debug/"
"Fanned_Single bar reduced model.ply");
// vcg::Matrix44d R;
// R.SetRotateDeg(45, {0, 0, 1});
// vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(mesh, R);
// mesh.updateEigenEdgeAndVertices();
FormFinder formFinder;
formFinder.runUnitTests();
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
fixedVertices[1] = std::unordered_set<DoFType>{2};
fixedVertices[2] = std::unordered_set<DoFType>{2};
fixedVertices[3] = std::unordered_set<DoFType>{1, 2};
fixedVertices[4] = std::unordered_set<DoFType>{2};
fixedVertices[5] = std::unordered_set<DoFType>{2};
fixedVertices[6] = std::unordered_set<DoFType>{0, 1, 2};
// Forced displacements
std::unordered_map<size_t, Eigen::Vector3d> nodalForcedDisplacements;
// nodalForcedDisplacements[10] = Eigen::Vector3d(0, 0.1, 0.1);
std::unordered_map<size_t, VectorType> nodalForcedNormal;
CoordType v14 = (mesh.vert[4].cP() - mesh.vert[0].cP()).Normalize();
CoordType v25 = (mesh.vert[5].cP() - mesh.vert[0].cP()).Normalize();
CoordType v36 = (mesh.vert[6].cP() - mesh.vert[0].cP()).Normalize();
const double forceMagnitude = 0.4;
std::unordered_map<VertexIndex, Vector6d> nodalForces{
// {4, Vector6d({0, 0, 0, v14[0], v14[1], 0}) * forceMagnitude},
// {1, Vector6d({0, 0, 0, -v14[0], -v14[1], 0}) * forceMagnitude},
// {5, Vector6d({0, 0, 0, v25[0], v25[1], 0}) * forceMagnitude},
// {2, Vector6d({0, 0, 0, -v25[0], -v25[1], 0}) * forceMagnitude},
// {6, Vector6d({0, 0, 0, v36[0], v36[1], 0}) * forceMagnitude},
{3, Vector6d({0, 0, 0, -v36[0], -v36[1], 0}) * forceMagnitude}};
// nodalForcedNormal[10] = v;
// nodalForcedNormal[0] = -v;
// fixedVertices[10] = {0, 1};
SimulationJob beamSimulationJob{std::make_shared<SimulationMesh>(mesh),
fixedVertices,
nodalForces,
{},
{}};
beamSimulationJob.mesh->setBeamMaterial(0.3, 1);
beamSimulationJob.mesh->setBeamCrossSection(
RectangularBeamDimensions{0.002, 0.002});
beamSimulationJob.registerForDrawing();
registerWorldAxes();
polyscope::show();
SimulationResults beamSimulationResults =
formFinder.executeSimulation(beamSimulationJob, true, true);
beamSimulationResults.simulationLabel = "Beam";
beamSimulationResults.registerForDrawing(beamSimulationJob);
polyscope::show();
// polyscope::show();
// SimulationResults beamSimulationResults =
// formFinder.executeSimulation(beamSimulationJob, true, true);
// beamSimulationResults.label = "Beam";
// beamSimulationResults.registerForDrawing(beamSimulationJob);
// polyscope::show();
return 0;
}

280
src/reducedmodeloptimizer.cpp
View File

@ -16,7 +16,8 @@ matplot::line_handle gPlotHandle;
std::vector<double> gObjectiveValueHistory;
Eigen::Vector4d g_initialX;
std::unordered_set<size_t> g_reducedPatternExludedEdges;
double g_initialParameters;
// double g_initialParameters;
Eigen::VectorXd g_initialParameters;
// struct OptimizationCallback {
// double operator()(const size_t &iterations, const Eigen::VectorXd &x,
@ -96,6 +97,7 @@ double g_initialParameters;
//};
double ReducedModelOptimizer::objective(long n, const double *x) {
std::cout.precision(17);
for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
std::cout << "x[" + std::to_string(parameterIndex) + "]="
@ -107,13 +109,14 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
// if (g_reducedPatternExludedEdges.contains(ei)) {
// continue;
// }
e.properties.E = g_initialParameters * x[0];
// e.properties.E = g_initialParameters * x[ei];
// e.properties.E = g_initialParameters(ei, 0) * x[ei];
// e.properties.G = g_initialParameters(1) * x[1];
// e.properties.A = g_initialParameters(0, 0) * x[0];
// e.properties.J = g_initialParameters(0, 1) * x[1];
// e.properties.I2 = g_initialParameters(0, 2) * x[2];
// e.properties.I3 = g_initialParameters(0, 3) * x[3];
e.properties.G = e.properties.E / (2 * (1 + 0.3));
e.properties.A = g_initialParameters(0) * x[0];
e.properties.J = g_initialParameters(1) * x[1];
e.properties.I2 = g_initialParameters(2) * x[2];
e.properties.I3 = g_initialParameters(3) * x[3];
// e.properties.G = e.properties.E / (2 * (1 + 0.3));
e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength;
e.torsionConstFactor = e.properties.G * e.properties.J / e.initialLength;
e.firstBendingConstFactor =
@ -145,10 +148,10 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
error += errorVector.norm();
}
// gObjectiveValueHistory.push_back(error);
// auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
// gObjectiveValueHistory.size());
// gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory);
gObjectiveValueHistory.push_back(error);
auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
gObjectiveValueHistory.size());
gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory);
return error;
}
@ -215,7 +218,7 @@ void ReducedModelOptimizer::computeMaps(
// reducedPattern.setLabel("ReducedPattern");
// Create opposite vertex map
m_fullPatternOppositeInterfaceViMap.clear();
for (size_t fanIndex = 0; fanIndex < fanSize / 2; fanIndex++) {
for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) {
const size_t vi0 = fullModelBaseTriangleInterfaceVi +
fanIndex * fullPatternInterfaceVertexOffset;
const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset;
@ -278,15 +281,19 @@ void ReducedModelOptimizer::computeMaps(
void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel,
FlatPattern &reducedModel) {
m_pReducedPatternElementalMesh =
if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) {
std::cerr << "A rectangular cross section is expected." << std::endl;
return;
}
m_pReducedPatternSimulationMesh =
std::make_shared<SimulationMesh>(reducedModel);
m_pReducedPatternElementalMesh->setBeamCrossSection(
m_pReducedPatternSimulationMesh->setBeamCrossSection(
CrossSectionType{0.002, 0.002});
m_pReducedPatternElementalMesh->setBeamMaterial(0.3, 1);
m_pFullModelElementalMesh = std::make_shared<SimulationMesh>(fullModel);
m_pFullModelElementalMesh->setBeamCrossSection(
m_pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1);
m_pFullModelSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
m_pFullModelSimulationMesh->setBeamCrossSection(
CrossSectionType{0.002, 0.002});
m_pFullModelElementalMesh->setBeamMaterial(0.3, 1);
m_pFullModelSimulationMesh->setBeamMaterial(0.3, 1);
}
ReducedModelOptimizer::ReducedModelOptimizer(
@ -312,7 +319,7 @@ void ReducedModelOptimizer::initialize(
}
void ReducedModelOptimizer::initializeStiffnesses() {
// g_initialParameters.resize(1, 2);
g_initialParameters.resize(4);
// Save save the beam stiffnesses
// for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) {
// Element &e = pReducedModelElementalMesh->elements[ei];
@ -323,18 +330,22 @@ void ReducedModelOptimizer::initializeStiffnesses() {
// e.firstBendingConstFactor *= stiffnessFactor;
// e.secondBendingConstFactor *= stiffnessFactor;
// }
g_initialParameters =
m_pReducedPatternElementalMesh->elements[0].properties.E;
// g_initialParameters =
// m_pReducedPatternSimulationMesh->elements[0].properties.E;
// for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) {
// g_initialParameters(ei, 0) =
// m_pReducedPatternSimulationMesh->elements[ei].properties.E;
// }
// g_initialParameters(1) =
// pReducedModelElementalMesh->elements[0].properties.G;
// g_initialParameters(0, 0) =
// pReducedModelElementalMesh->elements[0].properties.A;
// g_initialParameters(0, 1) =
// pReducedModelElementalMesh->elements[0].properties.J;
// g_initialParameters(0, 2) =
// pReducedModelElementalMesh->elements[0].properties.I2;
// g_initialParameters(0, 3) =
// pReducedModelElementalMesh->elements[0].properties.I3;
g_initialParameters(0) =
m_pReducedPatternSimulationMesh->elements[0].properties.A;
g_initialParameters(1) =
m_pReducedPatternSimulationMesh->elements[0].properties.J;
g_initialParameters(2) =
m_pReducedPatternSimulationMesh->elements[0].properties.I2;
g_initialParameters(3) =
m_pReducedPatternSimulationMesh->elements[0].properties.I3;
// }
}
@ -363,7 +374,7 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob(
fullModelNodalForcedRotation.second;
}
simulationJobOfReducedModel = SimulationJob{m_pReducedPatternElementalMesh,
simulationJobOfReducedModel = SimulationJob{m_pReducedPatternSimulationMesh,
reducedModelFixedVertices,
reducedModelNodalForces,
{},
@ -384,7 +395,7 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
optimalDisplacementsOfReducedModel.resize(0, 3);
optimalDisplacementsOfReducedModel.resize(
m_pReducedPatternElementalMesh->VN(), 3);
m_pReducedPatternSimulationMesh->VN(), 3);
optimalDisplacementsOfReducedModel.setZero(
optimalDisplacementsOfReducedModel.rows(),
optimalDisplacementsOfReducedModel.cols());
@ -407,10 +418,10 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob(
// polyscope::show();
gObjectiveValueHistory.clear();
fullPatternSimulationJob.mesh->savePly(
"Fanned_" + m_pFullModelElementalMesh->getLabel() + ".ply");
"Fanned_" + m_pFullModelSimulationMesh->getLabel() + ".ply");
SimulationResults fullModelResults =
simulator.executeSimulation(fullPatternSimulationJob, false, true);
fullModelResults.simulationLabel = "fullModel";
simulator.executeSimulation(fullPatternSimulationJob, false);
fullModelResults.label = "fullModel";
fullModelResults.registerForDrawing(fullPatternSimulationJob);
registerWorldAxes();
polyscope::show();
@ -420,24 +431,27 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob(
g_reducedPatternSimulationJob);
// gReducedPatternSimulationJob.registerForDrawing();
// polyscope::show();
fullModelResults.registerForDrawing(fullPatternSimulationJob);
polyscope::show();
double (*pObjectiveFunction)(long, const double *) = &objective;
const size_t n = 1;
// g_initialParameters.rows();
const size_t npt = ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
// const size_t n = m_pReducedPatternSimulationMesh->EN();
const size_t n = 4;
const size_t npt = 2 * n;
// ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2);
assert(npt <= 2 * n + 1 && "The choice of the number of interpolation "
"conditions is not recommended.");
// Set initial guess of solution
const double initialGuess = 1;
std::vector<double> x(n, initialGuess);
std::vector<double> x(n, 2);
if (!initialGuess.empty()) {
x = initialGuess;
} // {0.10000000000000 001, 2, 1.9999999971613847, 6.9560343643347764};
// {1, 5.9277};
// {0.0001, 2, 2.000000005047502, 1.3055270196964464};
// {initialGuess(0), initialGuess(1), initialGuess(2),
// initialGuess(3)};
const double xMin = 1;
const double xMax = 100;
const double xMin = 0.05;
const double xMax = 10;
// assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) {
// return d >= xMax || d <= xMin;
// }));
@ -474,7 +488,7 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob(
std::cout << "Final objective value:" << error << std::endl;
reducedModelOptimizedResults.simulationLabel = "reducedModel";
reducedModelOptimizedResults.label = "reducedModel";
reducedModelOptimizedResults.registerForDrawing(
g_reducedPatternSimulationJob);
polyscope::show();
@ -486,65 +500,79 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob(
return eigenX;
}
void ReducedModelOptimizer::setInitialGuess(std::vector<double> v) {
initialGuess = v;
}
std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
const std::shared_ptr<SimulationMesh> &pMesh) {
std::vector<SimulationJob> scenarios;
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
std::unordered_map<VertexIndex, Vector6d> nodalForces;
// NOTE: Assuming that the first interface node lays on the y axis
const double forceMagnitude = 1;
// Assuming the patterns lays on the x-y plane
const CoordType patternPlaneNormal(0, 0, 1);
// Make the first interface node lay on the x axis
const size_t fullPatternFirstInterfaceNodeIndex =
m_fullPatternOppositeInterfaceViMap.begin()->second;
CoordType fullPatternFirstInterfaceNodePosition =
m_pFullModelElementalMesh->vert[fullPatternFirstInterfaceNodeIndex].cP();
const vcg::Matrix33d R = vcg::RotationMatrix(
fullPatternFirstInterfaceNodePosition,
CoordType(fullPatternFirstInterfaceNodePosition.Norm(), 0, 0), false);
std::for_each(m_pFullModelElementalMesh->vert.begin(),
m_pFullModelElementalMesh->vert.end(), [&](auto &v) {
v.P() = R * v.P();
v.N() = R * v.N();
});
std::for_each(m_pReducedPatternElementalMesh->vert.begin(),
m_pReducedPatternElementalMesh->vert.end(), [&](auto &v) {
v.P() = R * v.P();
v.N() = R * v.N();
});
m_pFullModelElementalMesh->updateEigenEdgeAndVertices();
m_pReducedPatternElementalMesh->updateEigenEdgeAndVertices();
// //// Axial
// for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
// CoordType forceDirection =
// (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// .Normalize();
// nodalForces[viPair.first] = Vector6d({forceDirection[0],
// forceDirection[1],
// forceDirection[2], 0, 0, 0}) *
// forceMagnitude * 10;
// fixedVertices[viPair.second] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// const size_t fullPatternFirstInterfaceNodeIndex =
// m_fullPatternOppositeInterfaceViMap.begin()->second;
// CoordType fullPatternFirstInterfaceNodePosition =
// m_pFullModelSimulationMesh->vert[fullPatternFirstInterfaceNodeIndex].cP();
// CoordType centerOfMass(0, 0, 0);
// for (size_t vi = 0; vi < pMesh->VN(); vi++) {
// centerOfMass = centerOfMass + pMesh->vert[vi].P();
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
// centerOfMass /= pMesh->VN();
// vcg::tri::UpdatePosition<SimulationMesh>::Translate(
// *m_pFullModelSimulationMesh, -centerOfMass);
// vcg::tri::UpdatePosition<SimulationMesh>::Translate(
// *m_pReducedPatternSimulationMesh, centerOfMass);
// //// In-plane Bending
// fixedVertices.clear();
// nodalForces.clear();
// for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
// CoordType v =
// (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// .Normalize();
// CoordType forceDirection = (v ^ patternPlaneNormal).Normalize();
// nodalForces[viPair.first] = Vector6d({forceDirection[0],
// forceDirection[1],
// forceDirection[2], 0, 0, 0}) *
// 0.40 * forceMagnitude;
// fixedVertices[viPair.second] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
// const vcg::Matrix33d R = vcg::RotationMatrix(
// fullPatternFirstInterfaceNodePosition,
// CoordType(fullPatternFirstInterfaceNodePosition.Norm(), 0, 0), false);
// std::for_each(m_pFullModelSimulationMesh->vert.begin(),
// m_pFullModelSimulationMesh->vert.end(), [&](auto &v) {
// v.P() = R * v.P();
// v.N() = R * v.N();
// });
// std::for_each(m_pReducedPatternSimulationMesh->vert.begin(),
// m_pReducedPatternSimulationMesh->vert.end(), [&](auto &v) {
// v.P() = R * v.P();
// v.N() = R * v.N();
// });
// m_pFullModelSimulationMesh->updateEigenEdgeAndVertices();
// m_pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
//// Axial
for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
CoordType forceDirection =
(pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
.Normalize();
nodalForces[viPair.first] = Vector6d({forceDirection[0], forceDirection[1],
forceDirection[2], 0, 0, 0}) *
forceMagnitude * 10;
fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
//// In-plane Bending
fixedVertices.clear();
nodalForces.clear();
for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
CoordType v =
(pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
.Normalize();
CoordType forceDirection = (v ^ patternPlaneNormal).Normalize();
nodalForces[viPair.first] = Vector6d({forceDirection[0], forceDirection[1],
forceDirection[2], 0, 0, 0}) *
0.40 * forceMagnitude;
fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
// //// Torsion
// fixedVertices.clear();
@ -569,38 +597,38 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces});
// // Out - of - plane bending.Pull towards Z
// fixedVertices.clear();
// nodalForces.clear();
// for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
// nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) *
// 1; fixedVertices[viPair.second] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
//// Out - of - plane bending.Pull towards Z
fixedVertices.clear();
nodalForces.clear();
for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * 1;
fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
// //// Double using moments
// fixedVertices.clear();
// nodalForces.clear();
// for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
// viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) {
// const auto viPair = *viPairIt;
// if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
// fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
// fixedVertices[viPair.second] = std::unordered_set<DoFType>{1, 2};
// } else {
// fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
// fixedVertices[viPair.second] = std::unordered_set<DoFType>{2};
// }
// CoordType v =
// (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// .Normalize();
// nodalForces[viPair.first] =
// Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude * 0.1;
// nodalForces[viPair.second] =
// Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 0.1;
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
//// Double using moments
fixedVertices.clear();
nodalForces.clear();
for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) {
const auto viPair = *viPairIt;
if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 2};
} else {
fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
fixedVertices[viPair.second] = std::unordered_set<DoFType>{2};
}
CoordType v =
(pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
.Normalize();
nodalForces[viPair.first] =
Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude * 0.1;
nodalForces[viPair.second] =
Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 0.1;
}
scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
//// Saddle
fixedVertices.clear();
@ -633,9 +661,11 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
Eigen::VectorXd ReducedModelOptimizer::optimize() {
std::vector<SimulationJob> simulationJobs =
createScenarios(m_pFullModelElementalMesh);
createScenarios(m_pFullModelSimulationMesh);
std::vector<Eigen::VectorXd> results;
for (const SimulationJob &job : simulationJobs) {
for (int i = 0; i < simulationJobs.size(); i++) {
std::cout << "Running simulation job:" << i << std::endl;
const SimulationJob &job = simulationJobs[i];
polyscope::removeAllStructures();
auto result = optimizeForSimulationJob(job);
results.push_back(result);

7
src/reducedmodeloptimizer.hpp
View File

@ -11,14 +11,15 @@ using FullPatternVertexIndex = VertexIndex;
using ReducedPatternVertexIndex = VertexIndex;
class ReducedModelOptimizer {
std::shared_ptr<SimulationMesh> m_pReducedPatternElementalMesh;
std::shared_ptr<SimulationMesh> m_pFullModelElementalMesh;
std::shared_ptr<SimulationMesh> m_pReducedPatternSimulationMesh;
std::shared_ptr<SimulationMesh> m_pFullModelSimulationMesh;
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
m_fullToReducedInterfaceViMap;
std::unordered_map<FullPatternVertexIndex, FullPatternVertexIndex>
m_fullPatternOppositeInterfaceViMap;
std::unordered_map<size_t, size_t> nodeToSlot;
std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
std::vector<double> initialGuess;
public:
Eigen::VectorXd optimize();
@ -35,6 +36,8 @@ public:
void initialize(FlatPattern &fullPattern, FlatPattern &reducedPatterm,
const std::unordered_set<size_t> &reducedModelExcludedEges);
void setInitialGuess(std::vector<double> v);
private:
void computeDesiredReducedModelDisplacements(
const SimulationResults &fullModelResults,