From 9884cd175f85d3cb634ba62ceb726b5221759381 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Mon, 15 Mar 2021 19:56:14 +0200 Subject: [PATCH] The linear simulation model is used for optimizing the reduced model. --- src/linearsimulationmodel.cpp | 21 + src/linearsimulationmodel.hpp | 215 +++++ src/main.cpp | 12 +- src/reducedmodeloptimizer.cpp | 1680 ++++++++++++++++----------------- src/reducedmodeloptimizer.hpp | 78 +- 5 files changed, 1118 insertions(+), 888 deletions(-) create mode 100644 src/linearsimulationmodel.cpp create mode 100644 src/linearsimulationmodel.hpp diff --git a/src/linearsimulationmodel.cpp b/src/linearsimulationmodel.cpp new file mode 100644 index 0000000..6e663c6 --- /dev/null +++ b/src/linearsimulationmodel.cpp @@ -0,0 +1,21 @@ +#include "linearsimulationmodel.hpp" +#include +#include +//#include +#include +#include +#include + + + + + + + + + + + + + + diff --git a/src/linearsimulationmodel.hpp b/src/linearsimulationmodel.hpp new file mode 100644 index 0000000..38bb2d0 --- /dev/null +++ b/src/linearsimulationmodel.hpp @@ -0,0 +1,215 @@ +#ifndef LINEARSIMULATIONMODEL_HPP +#define LINEARSIMULATIONMODEL_HPP + +//#include "beam.hpp" +#include "simulationresult.hpp" +#include "threed_beam_fea.h" +#include +#include + +// 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 nodalForces; +// std::vector beamDimensions; +// std::vector beamMaterial; +//}; + +// struct SimulationResults { +// std::vector 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 { +public: + LinearSimulationModel(){ + + } + static std::vector + getFeaElements(const std::shared_ptr &job) { + const int numberOfEdges = job->pMesh->EN(); + std::vector elements(numberOfEdges); + for (int edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) { + const SimulationMesh::CoordType &evn0 = + job->pMesh->edge[edgeIndex].cV(0)->cN(); + const SimulationMesh::CoordType &evn1 = + job->pMesh->edge[edgeIndex].cV(1)->cN(); + const std::vector nAverageVector{(evn0[0] + evn1[0]) / 2, + (evn0[1] + evn1[1]) / 2, + (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); + 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); + } + + return elements; + } + static std::vector + getFeaNodes(const std::shared_ptr &job) { + const int numberOfNodes = job->pMesh->VN(); + std::vector feaNodes(numberOfNodes); + for (int vi = 0; vi < numberOfNodes; vi++) { + const CoordType &p = job->pMesh->vert[vi].cP(); + feaNodes[vi] = fea::Node(p[0], p[1], p[2]); + } + + return feaNodes; + } + static std::vector + getFeaFixedNodes(const std::shared_ptr &job) { + std::vector boundaryConditions; + boundaryConditions.reserve(job->constrainedVertices.size() * 6); + for (auto fixedVertex : job->constrainedVertices) { + const int vertexIndex = fixedVertex.first; + for (int dofIndex : fixedVertex.second) { + boundaryConditions.emplace_back( + fea::BC(vertexIndex, static_cast(dofIndex), 0)); + } + } + + return boundaryConditions; + } + static std::vector + getFeaNodalForces(const std::shared_ptr &job) { + std::vector nodalForces; + nodalForces.reserve(job->nodalExternalForces.size() * 6); + + for (auto nodalForce : job->nodalExternalForces) { + for (int dofIndex = 0; dofIndex < 6; dofIndex++) { + if (nodalForce.second[dofIndex] == 0) { + continue; + } + nodalForces.emplace_back( + fea::Force(nodalForce.first, dofIndex, nodalForce.second[dofIndex])); + } + } + + return nodalForces; + } + static SimulationResults getResults(const fea::Summary &feaSummary) { + SimulationResults results; + + results.executionTime = feaSummary.total_time_in_ms * 1000; + // displacements + results.displacements.resize(feaSummary.num_nodes); + for (int vi = 0; vi < feaSummary.num_nodes; vi++) { + results.displacements[vi] = Vector6d(feaSummary.nodal_displacements[vi]); + } + + // // Convert forces + // // Convert to vector of eigen matrices of the form force component-> per + // // Edge + // const int numDof = 6; + // const size_t numberOfEdges = feaSummary.element_forces.size(); + // edgeForces = + // std::vector(numDof, Eigen::VectorXd(2 * + // numberOfEdges)); + // for (gsl::index edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) { + // for (gsl::index forceComponentIndex = 0; forceComponentIndex < numDof; + // forceComponentIndex++) { + // (edgeForces[forceComponentIndex])(2 * edgeIndex) = + // feaSummary.element_forces[edgeIndex][forceComponentIndex]; + // (edgeForces[forceComponentIndex])(2 * edgeIndex + 1) = + // feaSummary.element_forces[edgeIndex][numDof + + // forceComponentIndex]; + // } + // } + return results; + } + + SimulationResults + executeSimulation(const std::shared_ptr &simulationJob) { + assert(simulationJob->pMesh->VN() != 0); + fea::Job job(getFeaNodes(simulationJob), getFeaElements(simulationJob)); + // printInfo(job); + + // create the default options + fea::Options opts; + opts.save_elemental_forces = false; + opts.save_nodal_displacements = false; + opts.save_nodal_forces = false; + opts.save_report = false; + opts.save_tie_forces = false; + // if (!elementalForcesOutputFilepath.empty()) { + // opts.save_elemental_forces = true; + // opts.elemental_forces_filename = elementalForcesOutputFilepath; + // } + // if (!nodalDisplacementsOutputFilepath.empty()) { + // opts.save_nodal_displacements = true; + // opts.nodal_displacements_filename = nodalDisplacementsOutputFilepath; + // } + + // have the program output status updates + opts.verbose = false; + + // form an empty vector of ties + std::vector ties; + + // also create an empty list of equations + std::vector equations; + auto fixedVertices = getFeaFixedNodes(simulationJob); + auto nodalForces = getFeaNodalForces(simulationJob); + fea::Summary feaResults = + fea::solve(job, fixedVertices, nodalForces, ties, equations, opts); + + SimulationResults results = getResults(feaResults); + results.job = simulationJob; + return results; + } + // SimulationResults getResults() const; + // void setResultsNodalDisplacementCSVFilepath(const std::string + // &outputPath); void setResultsElementalForcesCSVFilepath(const std::string + // &outputPath); + +private: + // std::string nodalDisplacementsOutputFilepath{"nodal_displacement.csv"}; + // std::string elementalForcesOutputFilepath{"elemental_forces.csv"}; + // SimulationResults results; + + static void printInfo(const fea::Job &job) { + std::cout << "Details regarding the fea::Job:" << std::endl; + std::cout << "Nodes:" << std::endl; + for (fea::Node n : job.nodes) { + std::cout << n << std::endl; + } + std::cout << "Elements:" << std::endl; + for (Eigen::Vector2i e : job.elems) { + std::cout << e << std::endl; + } + } +}; + +#endif // LINEARSIMULATIONMODEL_HPP diff --git a/src/main.cpp b/src/main.cpp index 068ac64..33606fe 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,7 +1,6 @@ #include "beamformfinder.hpp" #include "csvfile.hpp" #include "edgemesh.hpp" -#include "flatpattern.hpp" #include "reducedmodeloptimizer.hpp" #include "simulationhistoryplotter.hpp" #include "trianglepattterntopology.hpp" @@ -29,13 +28,13 @@ int main(int argc, char *argv[]) { // Populate the pattern pair to be optimized ////Full pattern const std::string filepath_fullPattern = argv[1]; - FlatPattern fullPattern(filepath_fullPattern); + PatternGeometry fullPattern(filepath_fullPattern); fullPattern.setLabel( std::filesystem::path(filepath_fullPattern).stem().string()); fullPattern.scale(0.03); ////Reduced pattern const std::string filepath_reducedPattern = argv[2]; - FlatPattern reducedPattern(filepath_reducedPattern); + PatternGeometry reducedPattern(filepath_reducedPattern); reducedPattern.setLabel( std::filesystem::path(filepath_reducedPattern).stem().string()); reducedPattern.scale(0.03); @@ -44,15 +43,16 @@ int main(int argc, char *argv[]) { ReducedModelOptimizer::xRange beamWidth{"B", 0.5, 1.5}; ReducedModelOptimizer::xRange beamDimensionsRatio{"bOverh", 0.7, 1.3}; ReducedModelOptimizer::xRange beamE{"E", 0.1, 1.9}; - ReducedModelOptimizer::xRange innerHexagonSize{"HS", 0.1, 0.9}; + ReducedModelOptimizer::xRange innerHexagonSize{"HexSize", 0.1, 0.8}; + ReducedModelOptimizer::xRange innerHexagonAngle{"HexAngle", -29.5, 29.5}; ReducedModelOptimizer::Settings settings_optimization; settings_optimization.xRanges = {beamWidth, beamDimensionsRatio, beamE, - innerHexagonSize}; + innerHexagonSize, innerHexagonAngle}; const bool input_numberOfFunctionCallsDefined = argc >= 4; settings_optimization.numberOfFunctionCalls = input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100; settings_optimization.normalizationStrategy = - ReducedModelOptimizer::Settings::NormalizationStrategy::NonNormalized; + ReducedModelOptimizer::Settings::NormalizationStrategy::Epsilon; settings_optimization.normalizationParameter = 0.0003; settings_optimization.solutionAccuracy = 0.01; diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index f397733..6cba086 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -1,38 +1,34 @@ #include "reducedmodeloptimizer.hpp" -#include "flatpattern.hpp" +#include "linearsimulationmodel.hpp" #include "simulationhistoryplotter.hpp" +#include "trianglepatterngeometry.hpp" #include "trianglepattterntopology.hpp" #include #include #include struct GlobalOptimizationVariables { - std::vector g_optimalReducedModelDisplacements; - std::vector> fullPatternDisplacements; - std::vector objectiveNormalizationValues; - std::vector> fullPatternSimulationJobs; - std::vector> reducedPatternSimulationJobs; - std::unordered_map - reducedToFullInterfaceViMap; - matplot::line_handle gPlotHandle; - std::vector gObjectiveValueHistory; - Eigen::Vector2d g_initialX; - std::unordered_set reducedPatternExludedEdges; - Eigen::VectorXd g_initialParameters; - std::vector - simulationScenarioIndices; - std::vector g_innerHexagonVectors{6, VectorType(0, 0, 0)}; - double g_innerHexagonInitialPos{0}; - bool optimizeInnerHexagonSize{false}; - std::vector firstOptimizationRoundResults; - int g_firstRoundIterationIndex{0}; - double minY{DBL_MAX}; - std::vector minX; - std::vector> failedSimulationsXRatio; - int numOfSimulationCrashes{false}; - int numberOfFunctionCalls{0}; - int numberOfOptimizationParameters{3}; - ReducedModelOptimizer::Settings optimizationSettings; + std::vector g_optimalReducedModelDisplacements; + std::vector> fullPatternDisplacements; + std::vector objectiveNormalizationValues; + std::vector> fullPatternSimulationJobs; + std::vector> reducedPatternSimulationJobs; + std::unordered_map + reducedToFullInterfaceViMap; + matplot::line_handle gPlotHandle; + std::vector gObjectiveValueHistory; + Eigen::VectorXd g_initialParameters; + std::vector + simulationScenarioIndices; + std::vector g_innerHexagonVectors{6, VectorType(0, 0, 0)}; + double innerHexagonInitialRotationAngle{30}; + double g_innerHexagonInitialPos{0}; + double minY{DBL_MAX}; + std::vector minX; + int numOfSimulationCrashes{false}; + int numberOfFunctionCalls{0}; + int numberOfOptimizationParameters{5}; + ReducedModelOptimizer::Settings optimizationSettings; } global; std::vector reducedPatternMaximumDisplacementSimulationJobs; @@ -43,14 +39,14 @@ double ReducedModelOptimizer::computeError( const std::unordered_map &reducedToFullInterfaceViMap, const double &normalizationFactor) { - const double rawError = - computeRawError(reducedPatternDisplacements, fullPatternDisplacements, - reducedToFullInterfaceViMap); - if (global.optimizationSettings.normalizationStrategy != - Settings::NormalizationStrategy::NonNormalized) { - return rawError / normalizationFactor; - } - return rawError; + const double rawError = + computeRawError(reducedPatternDisplacements, fullPatternDisplacements, + reducedToFullInterfaceViMap); + if (global.optimizationSettings.normalizationStrategy != + Settings::NormalizationStrategy::NonNormalized) { + return rawError / normalizationFactor; + } + return rawError; } double ReducedModelOptimizer::computeRawError( @@ -58,237 +54,244 @@ double ReducedModelOptimizer::computeRawError( const std::vector &fullPatternDisplacements, const std::unordered_map &reducedToFullInterfaceViMap) { - double error = 0; - for (const auto reducedFullViPair : reducedToFullInterfaceViMap) { - const VertexIndex reducedModelVi = reducedFullViPair.first; - Eigen::Vector3d reducedVertexDisplacement( - reducedPatternDisplacements[reducedModelVi][0], - reducedPatternDisplacements[reducedModelVi][1], - reducedPatternDisplacements[reducedModelVi][2]); - if (!std::isfinite(reducedVertexDisplacement[0]) || - !std::isfinite(reducedVertexDisplacement[1]) || - !std::isfinite(reducedVertexDisplacement[2])) { - std::cout << "Displacements are not finite" << std::endl; - std::terminate(); + double error = 0; + for (const auto reducedFullViPair : reducedToFullInterfaceViMap) { + const VertexIndex reducedModelVi = reducedFullViPair.first; + Eigen::Vector3d reducedVertexDisplacement( + reducedPatternDisplacements[reducedModelVi][0], + reducedPatternDisplacements[reducedModelVi][1], + reducedPatternDisplacements[reducedModelVi][2]); + if (!std::isfinite(reducedVertexDisplacement[0]) || + !std::isfinite(reducedVertexDisplacement[1]) || + !std::isfinite(reducedVertexDisplacement[2])) { + std::cout << "Displacements are not finite" << std::endl; + std::terminate(); + } + const VertexIndex fullModelVi = reducedFullViPair.second; + Eigen::Vector3d fullVertexDisplacement( + fullPatternDisplacements[fullModelVi][0], + fullPatternDisplacements[fullModelVi][1], + fullPatternDisplacements[fullModelVi][2]); + Eigen::Vector3d errorVector = + fullVertexDisplacement - reducedVertexDisplacement; + // error += errorVector.squaredNorm(); + error += errorVector.norm(); } - const VertexIndex fullModelVi = reducedFullViPair.second; - Eigen::Vector3d fullVertexDisplacement( - fullPatternDisplacements[fullModelVi][0], - fullPatternDisplacements[fullModelVi][1], - fullPatternDisplacements[fullModelVi][2]); - Eigen::Vector3d errorVector = - fullVertexDisplacement - reducedVertexDisplacement; - // error += errorVector.squaredNorm(); - error += errorVector.norm(); - } - return error; + return error; } void updateMesh(long n, const double *x) { - std::shared_ptr &pReducedPatternSimulationMesh = - global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]] - ->pMesh; - // const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0]; - // std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << " - // " - // << elem.firstBendingConstFactor << " " - // << elem.secondBendingConstFactor << std::endl; - for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { - Element &e = pReducedPatternSimulationMesh->elements[ei]; - // if (g_reducedPatternExludedEdges.contains(ei)) { - // continue; - // } - // e.properties.E = g_initialParameters * x[ei]; - // e.properties.E = g_initialParameters(0) * x[0]; - // e.properties.G = g_initialParameters(1) * x[1]; - e.setDimensions( - RectangularBeamDimensions(global.g_initialParameters(0) * x[0], - global.g_initialParameters(0) * x[0] / - (global.g_initialParameters(1) * x[1]))); - e.setMaterial(ElementMaterial(e.material.poissonsRatio, - global.g_initialParameters(2) * x[2])); - // 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 = - // 2 * e.properties.E * e.properties.I2 / e.initialLength; - // e.secondBendingConstFactor = - // 2 * e.properties.E * e.properties.I3 / e.initialLength; - } - - // std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << " - // " - // << elem.firstBendingConstFactor << " " - // << elem.secondBendingConstFactor << std::endl; - // const Element &e = pReducedPatternSimulationMesh->elements[0]; - // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " - // << e.firstBendingConstFactor << " " << - // e.secondBendingConstFactor - // << std::endl; - - if (global.optimizeInnerHexagonSize) { - assert(pReducedPatternSimulationMesh->EN() == 12); - for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN(); - vi += 2) { - pReducedPatternSimulationMesh->vert[vi].P() = - global.g_innerHexagonVectors[vi / 2] * x[n - 1]; + std::shared_ptr &pReducedPatternSimulationMesh = + global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]] + ->pMesh; + // const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0]; + // std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << " + // " + // << elem.firstBendingConstFactor << " " + // << elem.secondBendingConstFactor << std::endl; + for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { + Element &e = pReducedPatternSimulationMesh->elements[ei]; + // if (g_reducedPatternExludedEdges.contains(ei)) { + // continue; + // } + // e.properties.E = g_initialParameters * x[ei]; + // e.properties.E = g_initialParameters(0) * x[0]; + // e.properties.G = g_initialParameters(1) * x[1]; + e.setDimensions( + RectangularBeamDimensions(global.g_initialParameters(0) * x[0], + global.g_initialParameters(0) * x[0] / + (global.g_initialParameters(1) * x[1]))); + e.setMaterial(ElementMaterial(e.material.poissonsRatio, + global.g_initialParameters(2) * x[2])); + // 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 = + // 2 * e.properties.E * e.properties.I2 / e.initialLength; + // e.secondBendingConstFactor = + // 2 * e.properties.E * e.properties.I3 / e.initialLength; } - } - pReducedPatternSimulationMesh->reset(); - pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); + // std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << " + // " + // << elem.firstBendingConstFactor << " " + // << elem.secondBendingConstFactor << std::endl; + // const Element &e = pReducedPatternSimulationMesh->elements[0]; + // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " + // << e.firstBendingConstFactor << " " << + // e.secondBendingConstFactor + // << std::endl; + + assert(pReducedPatternSimulationMesh->EN() == 12); + auto R = vcg::RotationMatrix( + ReducedModelOptimizer::patternPlaneNormal, + vcg::math::ToRad(x[4] - global.innerHexagonInitialRotationAngle)); + // for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN(); + // vi += 2) { + for (int rotationCounter = 0; + rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) { + pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() = + R * global.g_innerHexagonVectors[rotationCounter] * x[3]; + } + + pReducedPatternSimulationMesh->reset(); +#ifdef POLYSCOPE_DEFINED + pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); +#endif } double ReducedModelOptimizer::objective(double b, double r, double E) { - std::vector x{b, r, E}; - return ReducedModelOptimizer::objective(x.size(), x.data()); + std::vector x{b, r, E}; + return ReducedModelOptimizer::objective(x.size(), x.data()); } double ReducedModelOptimizer::objective(double b, double h, double E, - double innerHexagonSize) { - std::vector x{b, h, E, innerHexagonSize}; - return ReducedModelOptimizer::objective(x.size(), x.data()); + double innerHexagonSize, + double innerHexagonRotationAngle) { + std::vector x{b, h, E, innerHexagonSize, innerHexagonRotationAngle}; + return ReducedModelOptimizer::objective(x.size(), x.data()); } double ReducedModelOptimizer::objective(long n, const double *x) { - // std::cout.precision(17); + // std::cout.precision(17); - // const Element &e = - // global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout << - // e.axialConstFactor << " " << e.torsionConstFactor << " " - // << e.firstBendingConstFactor << " " << - // e.secondBendingConstFactor - // << std::endl; - updateMesh(n, x); - // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " - // << e.firstBendingConstFactor << " " << - // e.secondBendingConstFactor - // << std::endl; + // const Element &e = + // global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout << + // e.axialConstFactor << " " << e.torsionConstFactor << " " + // << e.firstBendingConstFactor << " " << + // e.secondBendingConstFactor + // << std::endl; + updateMesh(n, x); + // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " + // << e.firstBendingConstFactor << " " << + // e.secondBendingConstFactor + // << std::endl; - // run simulations - double totalError = 0; - FormFinder simulator; - FormFinder::Settings simulationSettings; - // simulationSettings.shouldDraw = true; - for (const int simulationScenarioIndex : global.simulationScenarioIndices) { - // std::cout << "Executing " - // << simulationScenarioStrings[simulationScenarioIndex] - // << std::endl; - SimulationResults reducedModelResults = simulator.executeSimulation( - global.reducedPatternSimulationJobs[simulationScenarioIndex], - simulationSettings); - // std::string filename; - if (!reducedModelResults.converged) { - totalError += std::numeric_limits::max(); - global.numOfSimulationCrashes++; - bool beVerbose = false; // TODO: use an extern or data member for that - if (beVerbose) { - std::cout << "Failed simulation" << std::endl; - } - } else { - double simulationScenarioError = computeError( - reducedModelResults.displacements, - global.fullPatternDisplacements[simulationScenarioIndex], - global.reducedToFullInterfaceViMap, - global.objectiveNormalizationValues[simulationScenarioIndex]); + // run simulations + double totalError = 0; + LinearSimulationModel simulator; + for (const int simulationScenarioIndex : global.simulationScenarioIndices) { + SimulationResults reducedModelResults = simulator.executeSimulation( + global.reducedPatternSimulationJobs[simulationScenarioIndex]); + //#ifdef POLYSCOPE_DEFINED + // global.reducedPatternSimulationJobs[simulationScenarioIndex] + // ->pMesh->registerForDrawing(colors.reducedInitial); + // reducedModelResults.registerForDrawing(colors.reducedDeformed); + // polyscope::show(); + // reducedModelResults.unregister(); + //#endif + // std::string filename; + if (!reducedModelResults.converged) { + totalError += std::numeric_limits::max(); + global.numOfSimulationCrashes++; +#ifdef POLYSCOPE_DEFINED + std::cout << "Failed simulation" << std::endl; +#endif + } else { + double simulationScenarioError = computeError( + reducedModelResults.displacements, + global.fullPatternDisplacements[simulationScenarioIndex], + global.reducedToFullInterfaceViMap, + global.objectiveNormalizationValues[simulationScenarioIndex]); - // if (global.optimizationSettings.normalizationStrategy != - // NormalizationStrategy::Epsilon && - // simulationScenarioError > 1) { - // std::cout << "Simulation scenario " - // << - // simulationScenarioStrings[simulationScenarioIndex] - // << " results in an error bigger than one." << - // std::endl; - // for (size_t parameterIndex = 0; parameterIndex < n; - // parameterIndex++) { - // std::cout << "x[" + std::to_string(parameterIndex) + "]=" - // << x[parameterIndex] << std::endl; - // } - // } + // if (global.optimizationSettings.normalizationStrategy != + // NormalizationStrategy::Epsilon && + // simulationScenarioError > 1) { + // std::cout << "Simulation scenario " + // << + // simulationScenarioStrings[simulationScenarioIndex] + // << " results in an error bigger than one." << + // std::endl; + // for (size_t parameterIndex = 0; parameterIndex < n; + // parameterIndex++) { + // std::cout << "x[" + std::to_string(parameterIndex) + "]=" + // << x[parameterIndex] << std::endl; + // } + // } - //#ifdef POLYSCOPE_DEFINED - // ReducedModelOptimizer::visualizeResults( - // global.fullPatternSimulationJobs[simulationScenarioIndex], - // global.reducedPatternSimulationJobs[simulationScenarioIndex], - // global.reducedToFullInterfaceViMap, false); + //#ifdef POLYSCOPE_DEFINED + // ReducedModelOptimizer::visualizeResults( + // global.fullPatternSimulationJobs[simulationScenarioIndex], + // global.reducedPatternSimulationJobs[simulationScenarioIndex], + // global.reducedToFullInterfaceViMap, false); - // ReducedModelOptimizer::visualizeResults( - // global.fullPatternSimulationJobs[simulationScenarioIndex], - // std::make_shared( - // reducedPatternMaximumDisplacementSimulationJobs - // [simulationScenarioIndex]), - // global.reducedToFullInterfaceViMap, true); - // polyscope::removeAllStructures(); - //#endif // POLYSCOPE_DEFINED - totalError += simulationScenarioError; + // ReducedModelOptimizer::visualizeResults( + // global.fullPatternSimulationJobs[simulationScenarioIndex], + // std::make_shared( + // reducedPatternMaximumDisplacementSimulationJobs + // [simulationScenarioIndex]), + // global.reducedToFullInterfaceViMap, true); + // polyscope::removeAllStructures(); + //#endif // POLYSCOPE_DEFINED + totalError += simulationScenarioError; + } } - } - // std::cout << error << std::endl; - if (totalError < global.minY) { - global.minY = totalError; - global.minX.assign(x, x + n); - } - // if (++global.numberOfFunctionCalls %100== 0) { - // std::cout << "Number of function calls:" << global.numberOfFunctionCalls - // << std::endl; - //} + // std::cout << error << std::endl; + if (totalError < global.minY) { + global.minY = totalError; + global.minX.assign(x, x + n); + } +#ifdef POLYSCOPE_DEFINED + if (++global.numberOfFunctionCalls % 100 == 0) { + std::cout << "Number of function calls:" << global.numberOfFunctionCalls + << std::endl; + } +#endif + // compute error and return it + // global.gObjectiveValueHistory.push_back(error); + // auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), + // gObjectiveValueHistory.size()); + // std::vector colors(gObjectiveValueHistory.size(), 2); + // if (g_firstRoundIterationIndex != 0) { + // for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(), + // [](double &c) { c = 0.7; }); + // } + // gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors); + // SimulationResultsReporter::createPlot("Number of Steps", "Objective + // value", + // gObjectiveValueHistory); - // compute error and return it - // global.gObjectiveValueHistory.push_back(error); - // auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), - // gObjectiveValueHistory.size()); - // std::vector colors(gObjectiveValueHistory.size(), 2); - // if (g_firstRoundIterationIndex != 0) { - // for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(), - // [](double &c) { c = 0.7; }); - // } - // gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors); - // SimulationResultsReporter::createPlot("Number of Steps", "Objective - // value", - // gObjectiveValueHistory); - - return totalError; + return totalError; } void ReducedModelOptimizer::createSimulationMeshes( - FlatPattern &fullModel, FlatPattern &reducedModel, + PatternGeometry &fullModel, PatternGeometry &reducedModel, std::shared_ptr &pFullPatternSimulationMesh, std::shared_ptr &pReducedPatternSimulationMesh) { - if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) { - std::cerr << "Error: A rectangular cross section is expected." << std::endl; - terminate(); - } - // Full pattern - pFullPatternSimulationMesh = std::make_shared(fullModel); - pFullPatternSimulationMesh->setBeamCrossSection( - CrossSectionType{0.002, 0.002}); - pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); + if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) { + std::cerr << "Error: A rectangular cross section is expected." << std::endl; + terminate(); + } + // Full pattern + pFullPatternSimulationMesh = std::make_shared(fullModel); + pFullPatternSimulationMesh->setBeamCrossSection( + CrossSectionType{0.002, 0.002}); + pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); - // Reduced pattern - pReducedPatternSimulationMesh = - std::make_shared(reducedModel); - pReducedPatternSimulationMesh->setBeamCrossSection( - CrossSectionType{0.002, 0.002}); - pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); + // Reduced pattern + pReducedPatternSimulationMesh = + std::make_shared(reducedModel); + pReducedPatternSimulationMesh->setBeamCrossSection( + CrossSectionType{0.002, 0.002}); + pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); } -void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel, - FlatPattern &reducedModel) { - ReducedModelOptimizer::createSimulationMeshes( - fullModel, reducedModel, m_pFullPatternSimulationMesh, - m_pReducedPatternSimulationMesh); +void ReducedModelOptimizer::createSimulationMeshes( + PatternGeometry &fullModel, PatternGeometry &reducedModel) { + ReducedModelOptimizer::createSimulationMeshes( + fullModel, reducedModel, m_pFullPatternSimulationMesh, + m_pReducedPatternSimulationMesh); } void ReducedModelOptimizer::computeMaps( const std::unordered_set &reducedModelExcludedEdges, const std::unordered_map> &slotToNode, - FlatPattern &fullPattern, FlatPattern &reducedPattern, + PatternGeometry &fullPattern, PatternGeometry &reducedPattern, std::unordered_map &reducedToFullInterfaceViMap, std::unordered_map @@ -296,255 +299,255 @@ void ReducedModelOptimizer::computeMaps( std::unordered_map &fullPatternOppositeInterfaceViMap) { - // Compute the offset between the interface nodes - const size_t interfaceSlotIndex = 4; // bottom edge - assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() && - slotToNode.find(interfaceSlotIndex)->second.size() == 1); - // Assuming that in the bottom edge there is only one vertex which is also the - // interface - const size_t baseTriangleInterfaceVi = - *(slotToNode.find(interfaceSlotIndex)->second.begin()); + // Compute the offset between the interface nodes + const size_t interfaceSlotIndex = 4; // bottom edge + assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() && + slotToNode.find(interfaceSlotIndex)->second.size() == 1); + // Assuming that in the bottom edge there is only one vertex which is also the + // interface + const size_t baseTriangleInterfaceVi = + *(slotToNode.find(interfaceSlotIndex)->second.begin()); - vcg::tri::Allocator::PointerUpdater - pu_fullModel; - fullPattern.deleteDanglingVertices(pu_fullModel); - const size_t fullModelBaseTriangleInterfaceVi = - pu_fullModel.remap.empty() ? baseTriangleInterfaceVi - : pu_fullModel.remap[baseTriangleInterfaceVi]; - const size_t fullModelBaseTriangleVN = fullPattern.VN(); - fullPattern.createFan(); - const size_t duplicateVerticesPerFanPair = - fullModelBaseTriangleVN - fullPattern.VN() / 6; - const size_t fullPatternInterfaceVertexOffset = - fullModelBaseTriangleVN - duplicateVerticesPerFanPair; - // std::cout << "Dups in fan pair:" << duplicateVerticesPerFanPair << - // std::endl; + vcg::tri::Allocator::PointerUpdater< + PatternGeometry::VertexPointer> + pu_fullModel; + fullPattern.deleteDanglingVertices(pu_fullModel); + const size_t fullModelBaseTriangleInterfaceVi = + pu_fullModel.remap.empty() ? baseTriangleInterfaceVi + : pu_fullModel.remap[baseTriangleInterfaceVi]; + const size_t fullModelBaseTriangleVN = fullPattern.VN(); + fullPattern.createFan(); + const size_t duplicateVerticesPerFanPair = + fullModelBaseTriangleVN - fullPattern.VN() / 6; + const size_t fullPatternInterfaceVertexOffset = + fullModelBaseTriangleVN - duplicateVerticesPerFanPair; + // std::cout << "Dups in fan pair:" << duplicateVerticesPerFanPair << + // std::endl; - // Save excluded edges TODO:this changes the global object. Should this be a - // function parameter? - global.reducedPatternExludedEdges.clear(); - const size_t fanSize = 6; - const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN(); - for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { - for (const size_t ei : reducedModelExcludedEdges) { - global.reducedPatternExludedEdges.insert( - fanIndex * reducedBaseTriangleNumberOfEdges + ei); + // Save excluded edges TODO:this changes the global object. Should this be a + // function parameter? + // global.reducedPatternExludedEdges.clear(); + // const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN(); + // for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { + // for (const size_t ei : reducedModelExcludedEdges) { + // global.reducedPatternExludedEdges.insert( + // fanIndex * reducedBaseTriangleNumberOfEdges + ei); + // } + // } + + // Construct reduced->full and full->reduced interface vi map + reducedToFullInterfaceViMap.clear(); + vcg::tri::Allocator::PointerUpdater< + PatternGeometry::VertexPointer> + pu_reducedModel; + reducedPattern.deleteDanglingVertices(pu_reducedModel); + const size_t reducedModelBaseTriangleInterfaceVi = + pu_reducedModel.remap[baseTriangleInterfaceVi]; + const size_t reducedModelInterfaceVertexOffset = + reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/; + reducedPattern.createFan(); + for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { + reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex + + reducedModelBaseTriangleInterfaceVi] = + fullModelBaseTriangleInterfaceVi + + fanIndex * fullPatternInterfaceVertexOffset; } - } + fullToReducedInterfaceViMap.clear(); + constructInverseMap(reducedToFullInterfaceViMap, fullToReducedInterfaceViMap); - // Construct reduced->full and full->reduced interface vi map - reducedToFullInterfaceViMap.clear(); - vcg::tri::Allocator::PointerUpdater - pu_reducedModel; - reducedPattern.deleteDanglingVertices(pu_reducedModel); - const size_t reducedModelBaseTriangleInterfaceVi = - pu_reducedModel.remap[baseTriangleInterfaceVi]; - const size_t reducedModelInterfaceVertexOffset = - reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/; - reducedPattern.createFan(); - for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { - reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex + - reducedModelBaseTriangleInterfaceVi] = - fullModelBaseTriangleInterfaceVi + - fanIndex * fullPatternInterfaceVertexOffset; - } - fullToReducedInterfaceViMap.clear(); - constructInverseMap(reducedToFullInterfaceViMap, fullToReducedInterfaceViMap); + // Create opposite vertex map + fullPatternOppositeInterfaceViMap.clear(); + for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) { + const size_t vi0 = fullModelBaseTriangleInterfaceVi + + fanIndex * fullPatternInterfaceVertexOffset; + const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset; + assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN()); + fullPatternOppositeInterfaceViMap[vi0] = vi1; + } - // Create opposite vertex map - fullPatternOppositeInterfaceViMap.clear(); - for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) { - const size_t vi0 = fullModelBaseTriangleInterfaceVi + - fanIndex * fullPatternInterfaceVertexOffset; - const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset; - assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN()); - fullPatternOppositeInterfaceViMap[vi0] = vi1; - } - - const bool debugMapping = false; - if (debugMapping) { + const bool debugMapping = false; + if (debugMapping) { #if POLYSCOPE_DEFINED - reducedPattern.registerForDrawing(); - std::vector colors_reducedPatternExcludedEdges( - reducedPattern.EN(), glm::vec3(0, 0, 0)); - for (const size_t ei : global.reducedPatternExludedEdges) { - colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0); - } - const std::string label = reducedPattern.getLabel(); - polyscope::getCurveNetwork(label) - ->addEdgeColorQuantity("Excluded edges", - colors_reducedPatternExcludedEdges) - ->setEnabled(true); - polyscope::show(); + reducedPattern.registerForDrawing(); - std::vector nodeColorsOpposite(fullPattern.VN(), - glm::vec3(0, 0, 0)); - for (const std::pair oppositeVerts : - fullPatternOppositeInterfaceViMap) { - auto color = polyscope::getNextUniqueColor(); - nodeColorsOpposite[oppositeVerts.first] = color; - nodeColorsOpposite[oppositeVerts.second] = color; - } - fullPattern.registerForDrawing(); - polyscope::getCurveNetwork(fullPattern.getLabel()) - ->addNodeColorQuantity("oppositeMap", nodeColorsOpposite) - ->setEnabled(true); - polyscope::show(); + // std::vector colors_reducedPatternExcludedEdges( + // reducedPattern.EN(), glm::vec3(0, 0, 0)); + // for (const size_t ei : global.reducedPatternExludedEdges) { + // colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0); + // } + // const std::string label = reducedPattern.getLabel(); + // polyscope::getCurveNetwork(label) + // ->addEdgeColorQuantity("Excluded edges", + // colors_reducedPatternExcludedEdges) + // ->setEnabled(true); + // polyscope::show(); - std::vector nodeColorsReducedToFull_reduced(reducedPattern.VN(), - glm::vec3(0, 0, 0)); - std::vector nodeColorsReducedToFull_full(fullPattern.VN(), - glm::vec3(0, 0, 0)); - for (size_t vi = 0; vi < reducedPattern.VN(); vi++) { - if (global.reducedToFullInterfaceViMap.contains(vi)) { + std::vector nodeColorsOpposite(fullPattern.VN(), + glm::vec3(0, 0, 0)); + for (const std::pair oppositeVerts : + fullPatternOppositeInterfaceViMap) { + auto color = polyscope::getNextUniqueColor(); + nodeColorsOpposite[oppositeVerts.first] = color; + nodeColorsOpposite[oppositeVerts.second] = color; + } + fullPattern.registerForDrawing(); + polyscope::getCurveNetwork(fullPattern.getLabel()) + ->addNodeColorQuantity("oppositeMap", nodeColorsOpposite) + ->setEnabled(true); + polyscope::show(); - auto color = polyscope::getNextUniqueColor(); - nodeColorsReducedToFull_reduced[vi] = color; - nodeColorsReducedToFull_full[global.reducedToFullInterfaceViMap[vi]] = - color; - } - } - polyscope::getCurveNetwork(reducedPattern.getLabel()) - ->addNodeColorQuantity("reducedToFull_reduced", - nodeColorsReducedToFull_reduced) - ->setEnabled(true); - polyscope::getCurveNetwork(fullPattern.getLabel()) - ->addNodeColorQuantity("reducedToFull_full", - nodeColorsReducedToFull_full) - ->setEnabled(true); - polyscope::show(); + std::vector nodeColorsReducedToFull_reduced(reducedPattern.VN(), + glm::vec3(0, 0, 0)); + std::vector nodeColorsReducedToFull_full(fullPattern.VN(), + glm::vec3(0, 0, 0)); + for (size_t vi = 0; vi < reducedPattern.VN(); vi++) { + if (global.reducedToFullInterfaceViMap.contains(vi)) { + + auto color = polyscope::getNextUniqueColor(); + nodeColorsReducedToFull_reduced[vi] = color; + nodeColorsReducedToFull_full[global.reducedToFullInterfaceViMap[vi]] = + color; + } + } + polyscope::getCurveNetwork(reducedPattern.getLabel()) + ->addNodeColorQuantity("reducedToFull_reduced", + nodeColorsReducedToFull_reduced) + ->setEnabled(true); + polyscope::getCurveNetwork(fullPattern.getLabel()) + ->addNodeColorQuantity("reducedToFull_full", + nodeColorsReducedToFull_full) + ->setEnabled(true); + polyscope::show(); #endif - } + } } void ReducedModelOptimizer::computeMaps( - FlatPattern &fullPattern, FlatPattern &reducedPattern, + PatternGeometry &fullPattern, PatternGeometry &reducedPattern, const std::unordered_set &reducedModelExcludedEdges) { - ReducedModelOptimizer::computeMaps( - reducedModelExcludedEdges, slotToNode, fullPattern, reducedPattern, - global.reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap, - m_fullPatternOppositeInterfaceViMap); + ReducedModelOptimizer::computeMaps( + reducedModelExcludedEdges, slotToNode, fullPattern, reducedPattern, + global.reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap, + m_fullPatternOppositeInterfaceViMap); } ReducedModelOptimizer::ReducedModelOptimizer( const std::vector &numberOfNodesPerSlot) { - FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot); - FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode); + FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot); + FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode); } void ReducedModelOptimizer::initializePatterns( - FlatPattern &fullPattern, FlatPattern &reducedPattern, + PatternGeometry &fullPattern, PatternGeometry &reducedPattern, const std::unordered_set &reducedModelExcludedEdges) { - // fullPattern.setLabel("full_pattern_" + fullPattern.getLabel()); - // reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel()); - assert(fullPattern.VN() == reducedPattern.VN() && - fullPattern.EN() >= reducedPattern.EN()); + // fullPattern.setLabel("full_pattern_" + fullPattern.getLabel()); + // reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel()); + assert(fullPattern.VN() == reducedPattern.VN() && + fullPattern.EN() >= reducedPattern.EN()); #if POLYSCOPE_DEFINED - polyscope::removeAllStructures(); + polyscope::removeAllStructures(); #endif - // Create copies of the input models - FlatPattern copyFullPattern; - FlatPattern copyReducedPattern; - copyFullPattern.copy(fullPattern); - copyReducedPattern.copy(reducedPattern); - global.optimizeInnerHexagonSize = copyReducedPattern.EN() == 2; - if (global.optimizeInnerHexagonSize) { + // Create copies of the input models + PatternGeometry copyFullPattern; + PatternGeometry copyReducedPattern; + copyFullPattern.copy(fullPattern); + copyReducedPattern.copy(reducedPattern); + /* + * Here we create the vector that connects the central node with the bottom + * left node of the base triangle. During the optimization the vi%2==0 nodes + * move on these vectors. + * */ const double h = copyReducedPattern.getBaseTriangleHeight(); double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0)); VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h, 0); - - const int fanSize = 6; - const CoordType rotationAxis(0, 0, 1); - for (int rotationCounter = 0; rotationCounter < fanSize; - rotationCounter++) { - VectorType rotatedVector = - vcg::RotationMatrix(rotationAxis, - vcg::math::ToRad(rotationCounter * 60.0)) * - baseTriangle_leftBottomNode; - global.g_innerHexagonVectors[rotationCounter] = rotatedVector; + for (int rotationCounter = 0; rotationCounter < fanSize; rotationCounter++) { + VectorType rotatedVector = + vcg::RotationMatrix(patternPlaneNormal, + vcg::math::ToRad(rotationCounter * 60.0)) * + baseTriangle_leftBottomNode; + global.g_innerHexagonVectors[rotationCounter] = rotatedVector; } const double innerHexagonInitialPos_x = copyReducedPattern.vert[0].cP()[0] / global.g_innerHexagonVectors[0][0]; const double innerHexagonInitialPos_y = copyReducedPattern.vert[0].cP()[1] / global.g_innerHexagonVectors[0][1]; global.g_innerHexagonInitialPos = innerHexagonInitialPos_x; - } - computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges); - createSimulationMeshes(copyFullPattern, copyReducedPattern); - initializeOptimizationParameters(m_pReducedPatternSimulationMesh); + global.innerHexagonInitialRotationAngle = + 30; /* NOTE: Here I assume that the CW reduced pattern is given as input. + This is not very generic */ + computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges); + createSimulationMeshes(copyFullPattern, copyReducedPattern); + initializeOptimizationParameters(m_pReducedPatternSimulationMesh); } void ReducedModelOptimizer::initializeOptimizationParameters( const std::shared_ptr &mesh) { - global.numberOfOptimizationParameters = 3; - global.g_initialParameters.resize( - global.optimizeInnerHexagonSize ? ++global.numberOfOptimizationParameters - : global.numberOfOptimizationParameters); - // Save save the beam stiffnesses - // for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) { - // Element &e = pReducedModelElementalMesh->elements[ei]; - // if (g_reducedPatternExludedEdges.contains(ei)) { - // const double stiffnessFactor = 5; - // e.axialConstFactor *= stiffnessFactor; - // e.torsionConstFactor *= stiffnessFactor; - // e.firstBendingConstFactor *= stiffnessFactor; - // e.secondBendingConstFactor *= stiffnessFactor; - // } - const double initialB = std::sqrt(mesh->elements[0].A); - const double initialRatio = 1; - ; - global.g_initialParameters(0) = initialB; - global.g_initialParameters(1) = initialRatio; - global.g_initialParameters(2) = mesh->elements[0].material.youngsModulus; - if (global.optimizeInnerHexagonSize) { + global.numberOfOptimizationParameters = 5; + global.g_initialParameters.resize(global.numberOfOptimizationParameters); + // Save save the beam stiffnesses + // for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) { + // Element &e = pReducedModelElementalMesh->elements[ei]; + // if (g_reducedPatternExludedEdges.contains(ei)) { + // const double stiffnessFactor = 5; + // e.axialConstFactor *= stiffnessFactor; + // e.torsionConstFactor *= stiffnessFactor; + // e.firstBendingConstFactor *= stiffnessFactor; + // e.secondBendingConstFactor *= stiffnessFactor; + // } + const double initialB = std::sqrt(mesh->elements[0].A); + const double initialRatio = 1; + global.g_initialParameters(0) = initialB; + global.g_initialParameters(1) = initialRatio; + global.g_initialParameters(2) = mesh->elements[0].material.youngsModulus; global.g_initialParameters(3) = global.g_innerHexagonInitialPos; - } - // g_initialParameters = - // m_pReducedPatternSimulationMesh->elements[0].properties.E; - // for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) { - // } - // g_initialParameters(0) = mesh->elements[0].properties.E; - // g_initialParameters(1) = mesh->elements[0].properties.G; - // g_initialParameters(0) = mesh->elements[0].properties.A; - // g_initialParameters(1) = mesh->elements[0].properties.J; - // g_initialParameters(2) = mesh->elements[0].properties.I2; - // g_initialParameters(3) = mesh->elements[0].properties.I3; - // } + global.innerHexagonInitialRotationAngle = 30; + global.g_initialParameters(4) = global.innerHexagonInitialRotationAngle; + // g_initialParameters = + // m_pReducedPatternSimulationMesh->elements[0].properties.E; + // for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) { + // } + // g_initialParameters(0) = mesh->elements[0].properties.E; + // g_initialParameters(1) = mesh->elements[0].properties.G; + // g_initialParameters(0) = mesh->elements[0].properties.A; + // g_initialParameters(1) = mesh->elements[0].properties.J; + // g_initialParameters(2) = mesh->elements[0].properties.I2; + // g_initialParameters(3) = mesh->elements[0].properties.I3; + // } } void ReducedModelOptimizer::computeReducedModelSimulationJob( const SimulationJob &simulationJobOfFullModel, const std::unordered_map &simulationJobFullToReducedMap, SimulationJob &simulationJobOfReducedModel) { - assert(simulationJobOfReducedModel.pMesh->VN() != 0); - std::unordered_map> - reducedModelFixedVertices; - for (auto fullModelFixedVertex : - simulationJobOfFullModel.constrainedVertices) { - reducedModelFixedVertices[simulationJobFullToReducedMap.at( - fullModelFixedVertex.first)] = fullModelFixedVertex.second; - } + assert(simulationJobOfReducedModel.pMesh->VN() != 0); + std::unordered_map> + reducedModelFixedVertices; + for (auto fullModelFixedVertex : + simulationJobOfFullModel.constrainedVertices) { + reducedModelFixedVertices[simulationJobFullToReducedMap.at( + fullModelFixedVertex.first)] = fullModelFixedVertex.second; + } - std::unordered_map reducedModelNodalForces; - for (auto fullModelNodalForce : - simulationJobOfFullModel.nodalExternalForces) { - reducedModelNodalForces[simulationJobFullToReducedMap.at( - fullModelNodalForce.first)] = fullModelNodalForce.second; - } + std::unordered_map reducedModelNodalForces; + for (auto fullModelNodalForce : + simulationJobOfFullModel.nodalExternalForces) { + reducedModelNodalForces[simulationJobFullToReducedMap.at( + fullModelNodalForce.first)] = fullModelNodalForce.second; + } - // std::unordered_map - // reducedModelNodalForcedNormals; for (auto fullModelNodalForcedRotation : - // simulationJobOfFullModel.nodalForcedNormals) { - // reducedModelNodalForcedNormals[simulationJobFullToReducedMap.at( - // fullModelNodalForcedRotation.first)] = - // fullModelNodalForcedRotation.second; - // } - simulationJobOfReducedModel.constrainedVertices = reducedModelFixedVertices; - simulationJobOfReducedModel.nodalExternalForces = reducedModelNodalForces; - simulationJobOfReducedModel.label = simulationJobOfFullModel.getLabel(); - // simulationJobOfReducedModel.nodalForcedNormals = - // reducedModelNodalForcedNormals; + // std::unordered_map + // reducedModelNodalForcedNormals; for (auto fullModelNodalForcedRotation : + // simulationJobOfFullModel.nodalForcedNormals) { + // reducedModelNodalForcedNormals[simulationJobFullToReducedMap.at( + // fullModelNodalForcedRotation.first)] = + // fullModelNodalForcedRotation.second; + // } + simulationJobOfReducedModel.constrainedVertices = reducedModelFixedVertices; + simulationJobOfReducedModel.nodalExternalForces = reducedModelNodalForces; + simulationJobOfReducedModel.label = simulationJobOfFullModel.getLabel(); + // simulationJobOfReducedModel.nodalForcedNormals = + // reducedModelNodalForcedNormals; } #if POLYSCOPE_DEFINED @@ -556,59 +559,59 @@ void ReducedModelOptimizer::visualizeResults( const std::vector &simulationScenarios, const std::unordered_map &reducedToFullInterfaceViMap) { - FormFinder simulator; - std::shared_ptr pFullPatternSimulationMesh = - fullPatternSimulationJobs[0]->pMesh; - pFullPatternSimulationMesh->registerForDrawing(); - pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() + - "_undeformed.ply"); - reducedPatternSimulationJobs[0]->pMesh->savePly( - reducedPatternSimulationJobs[0]->pMesh->getLabel() + "_undeformed.ply"); - double totalError = 0; - for (const int simulationScenarioIndex : simulationScenarios) { - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationScenarioIndex]; - pFullPatternSimulationJob->registerForDrawing( - pFullPatternSimulationMesh->getLabel()); - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob); - fullModelResults.registerForDrawing(); - // fullModelResults.saveDeformedModel(); - const std::shared_ptr &pReducedPatternSimulationJob = - reducedPatternSimulationJobs[simulationScenarioIndex]; - SimulationResults reducedModelResults = - simulator.executeSimulation(pReducedPatternSimulationJob); - double normalizationFactor = 1; - if (global.optimizationSettings.normalizationStrategy != - Settings::NormalizationStrategy::NonNormalized) { - normalizationFactor = - global.objectiveNormalizationValues[simulationScenarioIndex]; + FormFinder simulator; + std::shared_ptr pFullPatternSimulationMesh = + fullPatternSimulationJobs[0]->pMesh; + pFullPatternSimulationMesh->registerForDrawing(); + pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() + + "_undeformed.ply"); + reducedPatternSimulationJobs[0]->pMesh->savePly( + reducedPatternSimulationJobs[0]->pMesh->getLabel() + "_undeformed.ply"); + double totalError = 0; + for (const int simulationScenarioIndex : simulationScenarios) { + const std::shared_ptr &pFullPatternSimulationJob = + fullPatternSimulationJobs[simulationScenarioIndex]; + pFullPatternSimulationJob->registerForDrawing( + pFullPatternSimulationMesh->getLabel()); + SimulationResults fullModelResults = + simulator.executeSimulation(pFullPatternSimulationJob); + fullModelResults.registerForDrawing(); + // fullModelResults.saveDeformedModel(); + const std::shared_ptr &pReducedPatternSimulationJob = + reducedPatternSimulationJobs[simulationScenarioIndex]; + SimulationResults reducedModelResults = + simulator.executeSimulation(pReducedPatternSimulationJob); + double normalizationFactor = 1; + if (global.optimizationSettings.normalizationStrategy != + Settings::NormalizationStrategy::NonNormalized) { + normalizationFactor = + global.objectiveNormalizationValues[simulationScenarioIndex]; + } + reducedModelResults.saveDeformedModel(); + fullModelResults.saveDeformedModel(); + double error = computeError( + reducedModelResults.displacements, fullModelResults.displacements, + reducedToFullInterfaceViMap, normalizationFactor); + std::cout << "Error of simulation scenario " + << simulationScenarioStrings[simulationScenarioIndex] << " is " + << error << std::endl; + totalError += error; + reducedModelResults.registerForDrawing(); + // firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing(); + // registerWorldAxes(); + const std::string screenshotFilename = + "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" + "Images/" + + pFullPatternSimulationMesh->getLabel() + "_" + + simulationScenarioStrings[simulationScenarioIndex]; + polyscope::show(); + polyscope::screenshot(screenshotFilename, false); + fullModelResults.unregister(); + reducedModelResults.unregister(); + // firstOptimizationRoundResults[simulationScenarioIndex].unregister(); } - reducedModelResults.saveDeformedModel(); - fullModelResults.saveDeformedModel(); - double error = computeError( - reducedModelResults.displacements, fullModelResults.displacements, - reducedToFullInterfaceViMap, normalizationFactor); - std::cout << "Error of simulation scenario " - << simulationScenarioStrings[simulationScenarioIndex] << " is " - << error << std::endl; - totalError += error; - reducedModelResults.registerForDrawing(); - // firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing(); - // registerWorldAxes(); - const std::string screenshotFilename = - "/home/iason/Coding/Projects/Approximating shapes with flat " - "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" - "Images/" + - pFullPatternSimulationMesh->getLabel() + "_" + - simulationScenarioStrings[simulationScenarioIndex]; - polyscope::show(); - polyscope::screenshot(screenshotFilename, false); - fullModelResults.unregister(); - reducedModelResults.unregister(); - // firstOptimizationRoundResults[simulationScenarioIndex].unregister(); - } - std::cout << "Total error:" << totalError << std::endl; + std::cout << "Total error:" << totalError << std::endl; } void ReducedModelOptimizer::registerResultsForDrawing( @@ -616,32 +619,32 @@ void ReducedModelOptimizer::registerResultsForDrawing( const std::shared_ptr &pReducedPatternSimulationJob, const std::unordered_map &reducedToFullInterfaceViMap) { - FormFinder simulator; - std::shared_ptr pFullPatternSimulationMesh = - pFullPatternSimulationJob->pMesh; - pFullPatternSimulationMesh->registerForDrawing(); - // pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() - // + - // "_undeformed.ply"); - // reducedPatternSimulationJobs[0]->pMesh->savePly( - // reducedPatternSimulationJobs[0]->pMesh->getLabel() + - // "_undeformed.ply"); - pFullPatternSimulationJob->registerForDrawing( - pFullPatternSimulationMesh->getLabel()); - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob); - fullModelResults.registerForDrawing(); - // fullModelResults.saveDeformedModel(); - SimulationResults reducedModelResults = - simulator.executeSimulation(pReducedPatternSimulationJob); - // reducedModelResults.saveDeformedModel(); - // fullModelResults.saveDeformedModel(); - double error = computeRawError( - reducedModelResults.displacements, fullModelResults.displacements, - reducedToFullInterfaceViMap/*, + FormFinder simulator; + std::shared_ptr pFullPatternSimulationMesh = + pFullPatternSimulationJob->pMesh; + pFullPatternSimulationMesh->registerForDrawing(); + // pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() + // + + // "_undeformed.ply"); + // reducedPatternSimulationJobs[0]->pMesh->savePly( + // reducedPatternSimulationJobs[0]->pMesh->getLabel() + + // "_undeformed.ply"); + pFullPatternSimulationJob->registerForDrawing( + pFullPatternSimulationMesh->getLabel()); + SimulationResults fullModelResults = + simulator.executeSimulation(pFullPatternSimulationJob); + fullModelResults.registerForDrawing(); + // fullModelResults.saveDeformedModel(); + SimulationResults reducedModelResults = + simulator.executeSimulation(pReducedPatternSimulationJob); + // reducedModelResults.saveDeformedModel(); + // fullModelResults.saveDeformedModel(); + double error = computeRawError( + reducedModelResults.displacements, fullModelResults.displacements, + reducedToFullInterfaceViMap/*, global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/); - std::cout << "Error is " << error << std::endl; - reducedModelResults.registerForDrawing(); + std::cout << "Error is " << error << std::endl; + reducedModelResults.registerForDrawing(); } #endif // POLYSCOPE_DEFINED @@ -650,396 +653,373 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( const std::unordered_map &displacementsReducedToFullMap, Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) { - assert(optimalDisplacementsOfReducedModel.rows() != 0 && - optimalDisplacementsOfReducedModel.cols() == 3); - optimalDisplacementsOfReducedModel.setZero( - optimalDisplacementsOfReducedModel.rows(), - optimalDisplacementsOfReducedModel.cols()); + assert(optimalDisplacementsOfReducedModel.rows() != 0 && + optimalDisplacementsOfReducedModel.cols() == 3); + optimalDisplacementsOfReducedModel.setZero( + optimalDisplacementsOfReducedModel.rows(), + optimalDisplacementsOfReducedModel.cols()); - for (auto reducedFullViPair : displacementsReducedToFullMap) { - const VertexIndex fullModelVi = reducedFullViPair.second; - const Vector6d fullModelViDisplacements = - fullModelResults.displacements[fullModelVi]; - optimalDisplacementsOfReducedModel.row(reducedFullViPair.first) = - Eigen::Vector3d(fullModelViDisplacements[0], - fullModelViDisplacements[1], - fullModelViDisplacements[2]); - } + for (auto reducedFullViPair : displacementsReducedToFullMap) { + const VertexIndex fullModelVi = reducedFullViPair.second; + const Vector6d fullModelViDisplacements = + fullModelResults.displacements[fullModelVi]; + optimalDisplacementsOfReducedModel.row(reducedFullViPair.first) = + Eigen::Vector3d(fullModelViDisplacements[0], + fullModelViDisplacements[1], + fullModelViDisplacements[2]); + } } ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(const Settings &settings) { - global.gObjectiveValueHistory.clear(); - dlib::matrix xMin(global.numberOfOptimizationParameters); - dlib::matrix xMax(global.numberOfOptimizationParameters); - for (int i = 0; i < global.numberOfOptimizationParameters; i++) { - xMin(i) = settings.xRanges[i].min; - xMax(i) = settings.xRanges[i].max; - } + global.gObjectiveValueHistory.clear(); + dlib::matrix xMin(global.numberOfOptimizationParameters); + dlib::matrix xMax(global.numberOfOptimizationParameters); + for (int i = 0; i < global.numberOfOptimizationParameters; i++) { + xMin(i) = settings.xRanges[i].min; + xMax(i) = settings.xRanges[i].max; + } - auto start = std::chrono::system_clock::now(); - dlib::function_evaluation result; - if (global.optimizeInnerHexagonSize) { - double (*objF)(double, double, double, double) = &objective; + auto start = std::chrono::system_clock::now(); + dlib::function_evaluation result; + double (*objF)(double, double, double, double, double) = &objective; result = dlib::find_min_global( objF, xMin, xMax, dlib::max_function_calls(settings.numberOfFunctionCalls), std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy); - } else { - double (*objF)(double, double, double) = &objective; - result = dlib::find_min_global( - objF, xMin, xMax, - dlib::max_function_calls(settings.numberOfFunctionCalls), - std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy); - } - auto end = std::chrono::system_clock::now(); - auto elapsed = - std::chrono::duration_cast(end - start); - Results results; - results.numberOfSimulationCrashes = global.numOfSimulationCrashes; - results.x = global.minX; - results.objectiveValue = global.minY; - if (global.minY != result.y) { - std::cerr << "Global min objective is not equal to result objective" - << std::endl; - } + auto end = std::chrono::system_clock::now(); + auto elapsed = + std::chrono::duration_cast(end - start); + Results results; + results.numberOfSimulationCrashes = global.numOfSimulationCrashes; + results.x = global.minX; + results.objectiveValue = global.minY; + if (global.minY != result.y) { + std::cerr << "Global min objective is not equal to result objective" + << std::endl; + } - // Compute raw objective value - if (global.optimizationSettings.normalizationStrategy != - Settings::NormalizationStrategy::NonNormalized) { - auto temp = global.optimizationSettings.normalizationStrategy; - global.optimizationSettings.normalizationStrategy = - Settings::NormalizationStrategy::NonNormalized; - if (global.optimizeInnerHexagonSize) { - results.rawObjectiveValue = objective(4, results.x.data()); + // Compute raw objective value + if (global.optimizationSettings.normalizationStrategy != + Settings::NormalizationStrategy::NonNormalized) { + auto temp = global.optimizationSettings.normalizationStrategy; + global.optimizationSettings.normalizationStrategy = + Settings::NormalizationStrategy::NonNormalized; + results.rawObjectiveValue = objective(results.x.size(), results.x.data()); + global.optimizationSettings.normalizationStrategy = temp; } else { - results.rawObjectiveValue = - objective(results.x[0], results.x[1], results.x[2]); + results.rawObjectiveValue = results.objectiveValue; } - global.optimizationSettings.normalizationStrategy = temp; - } else { - results.rawObjectiveValue = results.objectiveValue; - } + // Compute obj value per simulation scenario + updateMesh(results.x.size(), results.x.data()); + results.objectiveValuePerSimulationScenario.resize( + NumberOfSimulationScenarios); + FormFinder::Settings simulationSettings; + FormFinder simulator; + for (int simulationScenarioIndex = 0; + simulationScenarioIndex < NumberOfSimulationScenarios; + simulationScenarioIndex++) { + SimulationResults reducedModelResults = simulator.executeSimulation( + global.reducedPatternSimulationJobs[simulationScenarioIndex], + simulationSettings); + const double error = computeError( + reducedModelResults.displacements, + global.fullPatternDisplacements[simulationScenarioIndex], + global.reducedToFullInterfaceViMap, + global.objectiveNormalizationValues[simulationScenarioIndex]); - // Compute obj value per simulation scenario - updateMesh(results.x.size(), results.x.data()); - results.objectiveValuePerSimulationScenario.resize( - NumberOfSimulationScenarios); - FormFinder::Settings simulationSettings; - FormFinder simulator; - for (int simulationScenarioIndex = 0; - simulationScenarioIndex < NumberOfSimulationScenarios; - simulationScenarioIndex++) { - SimulationResults reducedModelResults = simulator.executeSimulation( - global.reducedPatternSimulationJobs[simulationScenarioIndex], - simulationSettings); - const double error = computeError( - reducedModelResults.displacements, - global.fullPatternDisplacements[simulationScenarioIndex], - global.reducedToFullInterfaceViMap, - global.objectiveNormalizationValues[simulationScenarioIndex]); + results.objectiveValuePerSimulationScenario[simulationScenarioIndex] = + error; + } + // if (result.y != + // std::accumulate(results.objectiveValuePerSimulationScenario.begin(), + // results.objectiveValuePerSimulationScenario.end(), 0)) + // { + // std::cerr + // << "Sum of per scenario objectives is not equal to result objective" + // << std::endl; + // } + results.time = elapsed.count() / 1000.0; + const bool printDebugInfo = false; + if (printDebugInfo) { + std::cout << "Finished optimizing." << endl; + // std::cout << "Solution x:" << endl; + // std::cout << result.x << endl; + std::cout << "Objective value:" << global.minY << endl; + } - results.objectiveValuePerSimulationScenario[simulationScenarioIndex] = - error; - } - // if (result.y != - // std::accumulate(results.objectiveValuePerSimulationScenario.begin(), - // results.objectiveValuePerSimulationScenario.end(), 0)) - // { - // std::cerr - // << "Sum of per scenario objectives is not equal to result objective" - // << std::endl; - // } - results.time = elapsed.count() / 1000.0; - const bool printDebugInfo = false; - if (printDebugInfo) { - std::cout << "Finished optimizing." << endl; - // std::cout << "Solution x:" << endl; - // std::cout << result.x << endl; - std::cout << "Objective value:" << global.minY << endl; - } - - return results; -} - -void ReducedModelOptimizer::setInitialGuess(std::vector v) { - initialGuess = v; + return results; } std::vector> ReducedModelOptimizer::createScenarios( const std::shared_ptr &pMesh) { - std::vector> scenarios; - scenarios.resize(SimulationScenario::NumberOfSimulationScenarios); - std::unordered_map> fixedVertices; - std::unordered_map nodalForces; - const double forceMagnitude = 1; - // Assuming the patterns lays on the x-y plane - const CoordType patternPlaneNormal(0, 0, 1); - // Assumes that the first interface node lays on the y axis + std::vector> scenarios; + scenarios.resize(SimulationScenario::NumberOfSimulationScenarios); + std::unordered_map> fixedVertices; + std::unordered_map nodalForces; + const double forceMagnitude = 1; - //// Axial - SimulationScenario scenarioName = SimulationScenario::Axial; - // NewMethod - for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); - viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { - if (viPairIt != m_fullPatternOppositeInterfaceViMap.begin()) { - CoordType forceDirection(1, 0, 0); - const auto viPair = *viPairIt; - nodalForces[viPair.first] = - Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, - 0, 0}) * - forceMagnitude * 2; - fixedVertices[viPair.second] = - std::unordered_set{0, 1, 2, 3, 4, 5}; + //// Axial + SimulationScenario scenarioName = SimulationScenario::Axial; + // NewMethod + for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); + viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { + if (viPairIt != m_fullPatternOppositeInterfaceViMap.begin()) { + CoordType forceDirection(1, 0, 0); + const auto viPair = *viPairIt; + nodalForces[viPair.first] = + Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, + 0, 0}) * + forceMagnitude * 2; + fixedVertices[viPair.second] = + std::unordered_set{0, 1, 2, 3, 4, 5}; + } } - } - // OldMethod - // 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{0, 1, 2, 3, 4, 5}; - // } - scenarios[scenarioName] = std::make_shared( - SimulationJob(pMesh, simulationScenarioStrings[scenarioName], - fixedVertices, nodalForces, {})); + // OldMethod + // 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{0, 1, 2, 3, 4, 5}; + // } + scenarios[scenarioName] = std::make_shared( + SimulationJob(pMesh, simulationScenarioStrings[scenarioName], + fixedVertices, nodalForces, {})); - //// Shear - scenarioName = SimulationScenario::Shear; - fixedVertices.clear(); - nodalForces.clear(); - // NewMethod - for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); - viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { - if (viPairIt != m_fullPatternOppositeInterfaceViMap.begin()) { - CoordType forceDirection(0, 1, 0); - const auto viPair = *viPairIt; - nodalForces[viPair.first] = - Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, - 0, 0}) * - forceMagnitude * 1; - fixedVertices[viPair.second] = - std::unordered_set{0, 1, 2, 3, 4, 5}; + //// Shear + scenarioName = SimulationScenario::Shear; + fixedVertices.clear(); + nodalForces.clear(); + // NewMethod + for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); + viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { + if (viPairIt != m_fullPatternOppositeInterfaceViMap.begin()) { + CoordType forceDirection(0, 1, 0); + const auto viPair = *viPairIt; + nodalForces[viPair.first] = + Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, + 0, 0}) * + forceMagnitude * 1; + fixedVertices[viPair.second] = + std::unordered_set{0, 1, 2, 3, 4, 5}; + } } - } - // OldMethod - // 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{0, 1, 2, 3, 4, 5}; - // } - scenarios[scenarioName] = std::make_shared( - SimulationJob(pMesh, simulationScenarioStrings[scenarioName], - fixedVertices, nodalForces, {})); + // OldMethod + // 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{0, 1, 2, 3, 4, 5}; + // } + scenarios[scenarioName] = std::make_shared( + SimulationJob(pMesh, simulationScenarioStrings[scenarioName], + fixedVertices, nodalForces, {})); - // //// Torsion - // fixedVertices.clear(); - // nodalForces.clear(); - // for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); - // viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { - // const auto &viPair = *viPairIt; - // if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { - // CoordType v = - // (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) - // .Normalize(); - // CoordType normalVDerivativeDir = (v ^ patternPlaneNormal).Normalize(); - // nodalForces[viPair.first] = Vector6d{ - // 0, 0, 0, normalVDerivativeDir[0], normalVDerivativeDir[1], 0}; - // fixedVertices[viPair.second] = - // std::unordered_set{0, 1, 2, 3, 4, 5}; - // fixedVertices[viPair.first] = std::unordered_set{0, 1, 2}; - // } else { - // fixedVertices[viPair.first] = std::unordered_set{0, 1, 2}; - // fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; - // } - // } - // scenarios.push_back({pMesh, fixedVertices, nodalForces}); + // //// Torsion + // fixedVertices.clear(); + // nodalForces.clear(); + // for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); + // viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { + // const auto &viPair = *viPairIt; + // if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { + // CoordType v = + // (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) + // .Normalize(); + // CoordType normalVDerivativeDir = (v ^ patternPlaneNormal).Normalize(); + // nodalForces[viPair.first] = Vector6d{ + // 0, 0, 0, normalVDerivativeDir[0], normalVDerivativeDir[1], 0}; + // fixedVertices[viPair.second] = + // std::unordered_set{0, 1, 2, 3, 4, 5}; + // fixedVertices[viPair.first] = std::unordered_set{0, 1, 2}; + // } else { + // fixedVertices[viPair.first] = std::unordered_set{0, 1, 2}; + // fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; + // } + // } + // scenarios.push_back({pMesh, fixedVertices, nodalForces}); - //// Bending - scenarioName = SimulationScenario::Bending; - 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{0, 1, 2, 3, 4, 5}; - } - scenarios[scenarioName] = std::make_shared( - SimulationJob(pMesh, simulationScenarioStrings[scenarioName], - fixedVertices, nodalForces, {})); - - //// Double using moments - scenarioName = SimulationScenario::Dome; - 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{0, 1, 2}; - fixedVertices[viPair.second] = std::unordered_set{0, 2}; - } else { - fixedVertices[viPair.first] = std::unordered_set{2}; - fixedVertices[viPair.second] = std::unordered_set{2}; + //// Bending + scenarioName = SimulationScenario::Bending; + 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{0, 1, 2, 3, 4, 5}; } - 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[scenarioName] = std::make_shared( - SimulationJob(pMesh, simulationScenarioStrings[scenarioName], - fixedVertices, nodalForces, {})); + scenarios[scenarioName] = std::make_shared( + SimulationJob(pMesh, simulationScenarioStrings[scenarioName], + fixedVertices, nodalForces, {})); - //// Saddle - scenarioName = SimulationScenario::Saddle; - fixedVertices.clear(); - nodalForces.clear(); - for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); - viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { - const auto &viPair = *viPairIt; - CoordType v = - (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) - .Normalize(); - if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { - nodalForces[viPair.first] = - Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.02 * forceMagnitude; - nodalForces[viPair.second] = - Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.02 * forceMagnitude; - } else { - fixedVertices[viPair.first] = std::unordered_set{2}; - fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; - - nodalForces[viPair.first] = - Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.01 * forceMagnitude; - nodalForces[viPair.second] = - Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.01 * forceMagnitude; + //// Double using moments + scenarioName = SimulationScenario::Dome; + 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{0, 1, 2}; + fixedVertices[viPair.second] = std::unordered_set{0, 2}; + } else { + fixedVertices[viPair.first] = std::unordered_set{2}; + fixedVertices[viPair.second] = std::unordered_set{2}; + } + CoordType v = + (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) ^ + CoordType(0, 0, -1).Normalize(); + nodalForces[viPair.first] = + Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude * 1; + nodalForces[viPair.second] = + Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 1; } - } - scenarios[scenarioName] = std::make_shared( - SimulationJob(pMesh, simulationScenarioStrings[scenarioName], - fixedVertices, nodalForces, {})); + scenarios[scenarioName] = std::make_shared( + SimulationJob(pMesh, simulationScenarioStrings[scenarioName], + fixedVertices, nodalForces, {})); - return scenarios; + //// Saddle + scenarioName = SimulationScenario::Saddle; + fixedVertices.clear(); + nodalForces.clear(); + for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); + viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { + const auto &viPair = *viPairIt; + CoordType v = + (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) ^ + CoordType(0, 0, -1).Normalize(); + if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { + nodalForces[viPair.first] = + Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.2 * forceMagnitude; + nodalForces[viPair.second] = + Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.2 * forceMagnitude; + } else { + fixedVertices[viPair.first] = std::unordered_set{2}; + fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; + + nodalForces[viPair.first] = + Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.1 * forceMagnitude; + nodalForces[viPair.second] = + Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.1 * forceMagnitude; + } + } + scenarios[scenarioName] = std::make_shared( + SimulationJob(pMesh, simulationScenarioStrings[scenarioName], + fixedVertices, nodalForces, {})); + + return scenarios; } void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() { - assert(global.optimizationSettings.shouldNormalizeObjectiveValue); - if (global.optimizationSettings.normalizationStrategy == - Settings::NormalizationStrategy::Epsilon) { + if (global.optimizationSettings.normalizationStrategy == + Settings::NormalizationStrategy::Epsilon) { - // Compute the sum of the displacement norms - std::vector fullPatternDisplacementNormSum( - NumberOfSimulationScenarios); - for (int simulationScenarioIndex : global.simulationScenarioIndices) { - double displacementNormSum = 0; - for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { - const Vector6d &vertexDisplacement = - global.fullPatternDisplacements[simulationScenarioIndex] - [interfaceViPair.second]; - displacementNormSum += vertexDisplacement.getTranslation().norm(); - } - fullPatternDisplacementNormSum[simulationScenarioIndex] = - displacementNormSum; - } - - for (int simulationScenarioIndex : global.simulationScenarioIndices) { - if (global.optimizationSettings.normalizationStrategy == - Settings::NormalizationStrategy::Epsilon) { - const double epsilon = - global.optimizationSettings.normalizationParameter; - if (epsilon > fullPatternDisplacementNormSum[simulationScenarioIndex]) { - // std::cout << "Epsilon used in " - // << - // simulationScenarioStrings[simulationScenarioIndex] - // << std::endl; + // Compute the sum of the displacement norms + std::vector fullPatternDisplacementNormSum( + NumberOfSimulationScenarios); + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + double displacementNormSum = 0; + for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { + const Vector6d &vertexDisplacement = + global.fullPatternDisplacements[simulationScenarioIndex] + [interfaceViPair.second]; + displacementNormSum += vertexDisplacement.getTranslation().norm(); + } + fullPatternDisplacementNormSum[simulationScenarioIndex] = + displacementNormSum; + } + + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + if (global.optimizationSettings.normalizationStrategy == + Settings::NormalizationStrategy::Epsilon) { + const double epsilon = + global.optimizationSettings.normalizationParameter; + if (epsilon > fullPatternDisplacementNormSum[simulationScenarioIndex]) { + // std::cout << "Epsilon used in " + // << + // simulationScenarioStrings[simulationScenarioIndex] + // << std::endl; + } + global.objectiveNormalizationValues[simulationScenarioIndex] = std::max( + fullPatternDisplacementNormSum[simulationScenarioIndex], epsilon); + // displacementNormSum; + } } - global.objectiveNormalizationValues[simulationScenarioIndex] = std::max( - fullPatternDisplacementNormSum[simulationScenarioIndex], epsilon); - // displacementNormSum; - } } - } } ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( const Settings &optimizationSettings, const std::vector &simulationScenarios) { - global.simulationScenarioIndices = simulationScenarios; - if (global.simulationScenarioIndices.empty()) { - global.simulationScenarioIndices = { - SimulationScenario::Axial, SimulationScenario::Shear, - SimulationScenario::Bending, SimulationScenario::Dome, - SimulationScenario::Saddle}; - } + global.simulationScenarioIndices = simulationScenarios; + if (global.simulationScenarioIndices.empty()) { + global.simulationScenarioIndices = { + SimulationScenario::Axial, SimulationScenario::Shear, + SimulationScenario::Bending, SimulationScenario::Dome, + SimulationScenario::Saddle}; + } - global.g_optimalReducedModelDisplacements.resize(NumberOfSimulationScenarios); - global.reducedPatternSimulationJobs.resize(NumberOfSimulationScenarios); - global.fullPatternDisplacements.resize(NumberOfSimulationScenarios); - global.objectiveNormalizationValues.resize(NumberOfSimulationScenarios); - global.g_firstRoundIterationIndex = 0; - global.minY = std::numeric_limits::max(); - global.numOfSimulationCrashes = 0; - global.numberOfFunctionCalls = 0; - global.optimizationSettings = optimizationSettings; - global.fullPatternSimulationJobs = - createScenarios(m_pFullPatternSimulationMesh); - reducedPatternMaximumDisplacementSimulationJobs.resize( - NumberOfSimulationScenarios); - // polyscope::removeAllStructures(); + global.g_optimalReducedModelDisplacements.resize(NumberOfSimulationScenarios); + global.reducedPatternSimulationJobs.resize(NumberOfSimulationScenarios); + global.fullPatternDisplacements.resize(NumberOfSimulationScenarios); + global.objectiveNormalizationValues.resize(NumberOfSimulationScenarios); + global.minY = std::numeric_limits::max(); + global.numOfSimulationCrashes = 0; + global.numberOfFunctionCalls = 0; + global.optimizationSettings = optimizationSettings; + global.fullPatternSimulationJobs = + createScenarios(m_pFullPatternSimulationMesh); + reducedPatternMaximumDisplacementSimulationJobs.resize( + NumberOfSimulationScenarios); + // polyscope::removeAllStructures(); - FormFinder::Settings simulationSettings; - // settings.shouldDraw = true; - for (int simulationScenarioIndex : global.simulationScenarioIndices) { - const std::shared_ptr &pFullPatternSimulationJob = - global.fullPatternSimulationJobs[simulationScenarioIndex]; - SimulationResults fullModelResults = simulator.executeSimulation( - pFullPatternSimulationJob, simulationSettings); - global.fullPatternDisplacements[simulationScenarioIndex] = - fullModelResults.displacements; - SimulationJob reducedPatternSimulationJob; - reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; - computeReducedModelSimulationJob(*pFullPatternSimulationJob, - m_fullToReducedInterfaceViMap, - reducedPatternSimulationJob); - global.reducedPatternSimulationJobs[simulationScenarioIndex] = - std::make_shared(reducedPatternSimulationJob); - } + FormFinder::Settings simulationSettings; + // settings.shouldDraw = true; + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + const std::shared_ptr &pFullPatternSimulationJob = + global.fullPatternSimulationJobs[simulationScenarioIndex]; + SimulationResults fullModelResults = simulator.executeSimulation( + pFullPatternSimulationJob, simulationSettings); + global.fullPatternDisplacements[simulationScenarioIndex] = + fullModelResults.displacements; + SimulationJob reducedPatternSimulationJob; + reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; + computeReducedModelSimulationJob(*pFullPatternSimulationJob, + m_fullToReducedInterfaceViMap, + reducedPatternSimulationJob); + global.reducedPatternSimulationJobs[simulationScenarioIndex] = + std::make_shared(reducedPatternSimulationJob); + } - if (global.optimizationSettings.normalizationStrategy != - Settings::NormalizationStrategy::NonNormalized) { - computeObjectiveValueNormalizationFactors(); - } - Results optResults = runOptimization(optimizationSettings); - for (int simulationScenarioIndex : global.simulationScenarioIndices) { - optResults.fullPatternSimulationJobs.push_back( - global.fullPatternSimulationJobs[simulationScenarioIndex]); - optResults.reducedPatternSimulationJobs.push_back( - global.reducedPatternSimulationJobs[simulationScenarioIndex]); - } + if (global.optimizationSettings.normalizationStrategy != + Settings::NormalizationStrategy::NonNormalized) { + computeObjectiveValueNormalizationFactors(); + } + Results optResults = runOptimization(optimizationSettings); + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + optResults.fullPatternSimulationJobs.push_back( + global.fullPatternSimulationJobs[simulationScenarioIndex]); + optResults.reducedPatternSimulationJobs.push_back( + global.reducedPatternSimulationJobs[simulationScenarioIndex]); + } #ifdef POLYSCOPE_DEFINED - optResults.draw(); + optResults.draw(); #endif // POLYSCOPE_DEFINED - return optResults; + return optResults; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index b9f759f..fd5514c 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -5,6 +5,7 @@ #include "csvfile.hpp" #include "edgemesh.hpp" #include "elementalmesh.hpp" +#include "linearsimulationmodel.hpp" #include "matplot/matplot.h" #include @@ -23,7 +24,6 @@ class ReducedModelOptimizer { m_fullPatternOppositeInterfaceViMap; std::unordered_map nodeToSlot; std::unordered_map> slotToNode; - std::vector initialGuess; #ifdef POLYSCOPE_DEFINED struct StaticColors { glm::vec3 fullInitial; @@ -41,6 +41,8 @@ class ReducedModelOptimizer { #endif // POLYSCOPE_DEFINED public: + inline static int fanSize{6}; + inline static VectorType patternPlaneNormal{0, 0, 1}; enum SimulationScenario { Axial, Shear, @@ -49,9 +51,8 @@ public: Saddle, NumberOfSimulationScenarios }; - - struct xRange { - std::string label; + struct xRange{ + std::string label; double min; double max; std::string toString() const { @@ -62,9 +63,14 @@ public: struct Results; struct Settings { - enum NormalizationStrategy { NonNormalized, Epsilon }; + enum NormalizationStrategy { + NonNormalized, + Epsilon, + MaxDisplacement, + EqualDisplacements + }; inline static vector normalizationStrategyStrings{ - "NonNormalized", "Epsilon"}; + "NonNormalized", "Epsilon", "MaxDsiplacement", "EqualDisplacements"}; std::vector xRanges; int numberOfFunctionCalls{100}; double solutionAccuracy{1e-2}; @@ -110,11 +116,8 @@ public: } os << numberOfFunctionCalls; os << solutionAccuracy; - if (normalizationStrategy == Epsilon) { - os << "Epsilon_" + std::to_string(normalizationParameter); - } else { - os << "NonNormalized"; - } + os << normalizationStrategyStrings[normalizationStrategy] + "_" + + std::to_string(normalizationParameter); } }; struct Results { @@ -183,7 +186,7 @@ public: const auto filepath = fileEntry.path(); if (filepath.extension() == ".json") { SimulationJob job; - job.load(filepath.string()); + job.load(filepath.string()); reducedPatternSimulationJobs.push_back( std::make_shared(job)); } @@ -194,12 +197,13 @@ public: void draw() const { initPolyscope(); FormFinder simulator; + LinearSimulationModel linearSimulator; assert(fullPatternSimulationJobs.size() == reducedPatternSimulationJobs.size()); fullPatternSimulationJobs[0]->pMesh->registerForDrawing( colors.fullInitial); reducedPatternSimulationJobs[0]->pMesh->registerForDrawing( - colors.reducedInitial); + colors.reducedInitial, false); const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); for (int simulationJobIndex = 0; @@ -212,14 +216,24 @@ public: SimulationResults fullModelResults = simulator.executeSimulation(pFullPatternSimulationJob); fullModelResults.registerForDrawing(colors.fullDeformed); + SimulationResults fullModelLinearResults = + linearSimulator.executeSimulation(pFullPatternSimulationJob); + fullModelLinearResults.setLabelPrefix("linear"); + fullModelLinearResults.registerForDrawing(colors.fullDeformed, false); // Drawing of reduced pattern results const std::shared_ptr &pReducedPatternSimulationJob = reducedPatternSimulationJobs[simulationJobIndex]; SimulationResults reducedModelResults = simulator.executeSimulation(pReducedPatternSimulationJob); reducedModelResults.registerForDrawing(colors.reducedDeformed); + SimulationResults reducedModelLinearResults = + linearSimulator.executeSimulation(pReducedPatternSimulationJob); + reducedModelLinearResults.setLabelPrefix("linear"); + reducedModelLinearResults.registerForDrawing(colors.reducedDeformed, + false); polyscope::options::programName = fullPatternSimulationJobs[0]->pMesh->getLabel(); + polyscope::view::resetCameraToHomeView(); polyscope::show(); // Save a screensh const std::string screenshotFilename = @@ -230,12 +244,14 @@ public: polyscope::screenshot(screenshotFilename, false); fullModelResults.unregister(); reducedModelResults.unregister(); + reducedModelLinearResults.unregister(); + fullModelLinearResults.unregister(); // double error = computeError( // reducedModelResults.displacements,fullModelResults.displacements, // ); // std::cout << "Error of simulation scenario " // << - // simulationScenarioStrings[simulationScenarioIndex] + // simula simulationScenarioStrings[simulationScenarioIndex] // << " is " // << error << std::endl; } @@ -333,31 +349,29 @@ public: getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); void initializePatterns( - FlatPattern &fullPattern, FlatPattern &reducedPatterm, + PatternGeometry &fullPattern, PatternGeometry &reducedPatterm, const std::unordered_set &reducedModelExcludedEges); - void setInitialGuess(std::vector v); - - static void runBeamOptimization(); - static void runSimulation(const std::string &filename, std::vector &x); - static double objective(double x0, double x1, double x2, double x3); + static double objective(double x0, double x1, double x2, double x3, + double innerHexagonRotationAngle); static double objective(double b, double r, double E); static std::vector> - createScenarios(const std::shared_ptr &pMesh, - const std::unordered_map - &fullPatternOppositeInterfaceViMap); + createScenarios(const std::shared_ptr &pMesh, + const std::unordered_map + &fullPatternOppositeInterfaceViMap); + static void createSimulationMeshes( - FlatPattern &fullModel, FlatPattern &reducedModel, + PatternGeometry &fullModel, PatternGeometry &reducedModel, std::shared_ptr &pFullPatternSimulationMesh, std::shared_ptr &pReducedPatternSimulationMesh); static void computeMaps( const std::unordered_set &reducedModelExcludedEdges, const std::unordered_map> &slotToNode, - FlatPattern &fullPattern, FlatPattern &reducedPattern, + PatternGeometry &fullPattern, PatternGeometry &reducedPattern, std::unordered_map &reducedToFullInterfaceViMap, std::unordered_map @@ -396,17 +410,17 @@ public: const double &normalizationFactor); private: - static void computeDesiredReducedModelDisplacements( - const SimulationResults &fullModelResults, - const std::unordered_map &displacementsReducedToFullMap, - Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); + static void computeDesiredReducedModelDisplacements( + const SimulationResults &fullModelResults, + const std::unordered_map &displacementsReducedToFullMap, + Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); static Results runOptimization(const Settings &settings); std::vector> createScenarios(const std::shared_ptr &pMesh); - void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern, + void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern, const std::unordered_set &reducedModelExcludedEges); - void createSimulationMeshes(FlatPattern &fullModel, - FlatPattern &reducedModel); + void createSimulationMeshes(PatternGeometry &fullModel, + PatternGeometry &reducedModel); static void initializeOptimizationParameters(const std::shared_ptr &mesh);