diff --git a/CMakeLists.txt b/CMakeLists.txt index 1229e3e..8ddb879 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -68,15 +68,6 @@ download_project(PROJ vcglib_devel ) file(GLOB EXT_SOURCES ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp) -##threed-beam-fea -download_project(PROJ threed-beam-fea - GIT_REPOSITORY https://github.com/IasonManolas/threed-beam-fea.git - GIT_TAG master - PREFIX ${EXTERNAL_DEPS_DIR} - ${UPDATE_DISCONNECTED_IF_AVAILABLE} -) -add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR}) - ##Eigen 3 NOTE: Eigen is required on the system the code is ran find_package(Eigen3 3.3 REQUIRED) @@ -94,7 +85,6 @@ target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20) target_include_directories(${PROJECT_NAME} PUBLIC ${vcglib_devel_SOURCE_DIR} - PUBLIC ${threed-beam-fea_SOURCE_DIR}/include PUBLIC ${MYSOURCES_SOURCE_DIR} PUBLIC ${MYSOURCES_SOURCE_DIR}/boost_graph ) @@ -102,7 +92,7 @@ target_include_directories(${PROJECT_NAME} target_link_directories(${PROJECT_NAME} PRIVATE ${MYSOURCES_SOURCE_DIR}/boost_graph/libs/) if(${USE_POLYSCOPE}) - target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen dlib::dlib ThreedBeamFEA MySources) + target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen dlib::dlib MySources) else() - target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen dlib::dlib ThreedBeamFEA MySources) + target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen dlib::dlib MySources) endif() diff --git a/src/linearsimulationmodel.cpp b/src/linearsimulationmodel.cpp deleted file mode 100644 index e69de29..0000000 diff --git a/src/linearsimulationmodel.hpp b/src/linearsimulationmodel.hpp deleted file mode 100644 index 38bb2d0..0000000 --- a/src/linearsimulationmodel.hpp +++ /dev/null @@ -1,215 +0,0 @@ -#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 a8eaf8e..3179038 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,4 +1,4 @@ -#include "beamformfinder.hpp" +#include "drmsimulationmodel.hpp" #include "csvfile.hpp" #include "edgemesh.hpp" #include "reducedmodeloptimizer.hpp" @@ -41,21 +41,21 @@ int main(int argc, char *argv[]) { reducedPattern.scale(0.03,interfaceNodeIndex); // Set the optization settings - ReducedModelOptimizer::xRange beamE{"E", 0.1, 1.9}; - ReducedModelOptimizer::xRange beamA{"A", 0.25, 2.25}; - ReducedModelOptimizer::xRange beamI2{"I2", std::sqrt(beamA.min)*std::pow(std::sqrt(beamA.min),3)/12, std::sqrt(beamA.max)*std::pow(std::sqrt(beamA.max),3)/12}; - ReducedModelOptimizer::xRange beamI3{"I3", std::sqrt(beamA.min)*std::pow(std::sqrt(beamA.min),3)/12, std::sqrt(beamA.max)*std::pow(std::sqrt(beamA.max),3)/12}; - ReducedModelOptimizer::xRange beamJ{"J", beamI2.min+beamI3.min, beamI2.max+beamI3.max}; - ReducedModelOptimizer::xRange innerHexagonSize{"HexSize", 0.1, 0.8}; - ReducedModelOptimizer::xRange innerHexagonAngle{"HexAngle", -29.5, 29.5}; - ReducedModelOptimizer::Settings settings_optimization; + ReducedModelOptimization::xRange beamE{"E", 0.001, 10000}; + ReducedModelOptimization::xRange beamA{"A", 0.001, 10000}; + ReducedModelOptimization::xRange beamI2{"I2", 0.001,10000}; + ReducedModelOptimization::xRange beamI3{"I3", 0.001,10000}; + ReducedModelOptimization::xRange beamJ{"J", 0.001,10000}; + ReducedModelOptimization::xRange innerHexagonSize{"HexSize", 0.1, 0.8}; + ReducedModelOptimization::xRange innerHexagonAngle{"HexAngle", -29.5, 29.5}; + ReducedModelOptimization::Settings settings_optimization; settings_optimization.xRanges = {beamE,beamA,beamJ,beamI2,beamI3, innerHexagonSize, innerHexagonAngle}; const bool input_numberOfFunctionCallsDefined = argc >= 4; settings_optimization.numberOfFunctionCalls = input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100; settings_optimization.normalizationStrategy = - ReducedModelOptimizer::Settings::NormalizationStrategy::Epsilon; + ReducedModelOptimization::Settings::NormalizationStrategy::Epsilon; settings_optimization.normalizationParameter = 0.0003; settings_optimization.solutionAccuracy = 0.01; @@ -66,8 +66,8 @@ int main(int argc, char *argv[]) { const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; assert(interfaceNodeIndex==numberOfNodesPerSlot[0]+numberOfNodesPerSlot[3]); ReducedModelOptimizer optimizer(numberOfNodesPerSlot); - optimizer.initializePatterns(fullPattern, reducedPattern, {}); - ReducedModelOptimizer::Results optimizationResults = + optimizer.initializePatterns(fullPattern, reducedPattern, {},settings_optimization.xRanges.size()); + ReducedModelOptimization::Results optimizationResults = optimizer.optimize(settings_optimization); // Export results diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index ee6815a..e0cfe5d 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -7,6 +7,8 @@ #include #include +using namespace ReducedModelOptimization; + struct GlobalOptimizationVariables { std::vector g_optimalReducedModelDisplacements; std::vector> fullPatternDisplacements; @@ -16,9 +18,9 @@ struct GlobalOptimizationVariables { std::unordered_map reducedToFullInterfaceViMap; matplot::line_handle gPlotHandle; - std::vector gObjectiveValueHistory; + std::vector objectiveValueHistory; Eigen::VectorXd initialParameters; - std::vector + std::vector simulationScenarioIndices; std::vector g_innerHexagonVectors{6, VectorType(0, 0, 0)}; double innerHexagonInitialRotationAngle{30}; @@ -28,7 +30,7 @@ struct GlobalOptimizationVariables { int numOfSimulationCrashes{false}; int numberOfFunctionCalls{0}; int numberOfOptimizationParameters{5}; - ReducedModelOptimizer::Settings optimizationSettings; + ReducedModelOptimization::Settings optimizationSettings; } global; std::vector reducedPatternMaximumDisplacementSimulationJobs; @@ -81,47 +83,6 @@ double ReducedModelOptimizer::computeRawError( return error; } -void updateMesh(long n, const double *x) { - std::shared_ptr &pReducedPatternSimulationMesh = - global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]] - ->pMesh; - - const double E=global.initialParameters(0)*x[0]; - const double A=global.initialParameters(1) * x[1]; - const double beamWidth=std::sqrt(A); - const double beamHeight=beamWidth; - -const double J=global.initialParameters(2) * x[2]; -const double I2=global.initialParameters(3) * x[3]; -const double I3=global.initialParameters(4) * x[4]; - for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { - Element &e = pReducedPatternSimulationMesh->elements[ei]; - e.setDimensions( - RectangularBeamDimensions(beamWidth, - beamHeight)); - e.setMaterial(ElementMaterial(e.material.poissonsRatio, - E)); - e.A = A; - e.J = J; - e.I2 = I2; - e.I3 = I3; - } - assert(pReducedPatternSimulationMesh->EN() == 12); - auto R = vcg::RotationMatrix( - ReducedModelOptimizer::patternPlaneNormal, - vcg::math::ToRad(x[6] - global.innerHexagonInitialRotationAngle)); - for (int rotationCounter = 0; - rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) { - pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() = - R * global.g_innerHexagonVectors[rotationCounter] * x[5]; - } - - 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()); @@ -134,8 +95,19 @@ double ReducedModelOptimizer::objective(double E,double A,double J,double I2,dou return ReducedModelOptimizer::objective(x.size(), x.data()); } +double ReducedModelOptimizer::objective(double E,double A, + double innerHexagonSize, + double innerHexagonRotationAngle) { + std::vector x{E,A, innerHexagonSize, innerHexagonRotationAngle}; + return ReducedModelOptimizer::objective(x.size(), x.data()); +} + double ReducedModelOptimizer::objective(long n, const double *x) { // std::cout.precision(17); +// for(int i=0;ipMesh->elements[0]; std::cout << @@ -219,19 +191,18 @@ double ReducedModelOptimizer::objective(long n, const double *x) { << 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.objectiveValueHistory.push_back(totalError); +// auto xPlot = matplot::linspace(0, global.objectiveValueHistory.size(), +// global.objectiveValueHistory.size()); +// std::vector colors(global.gObjectiveValueHistory.size(), 2); +// if (global.g_firstRoundIterationIndex != 0) { +// for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(), +// [](double &c) { c = 0.7; }); +// } +// global.gPlotHandle = matplot::scatter(xPlot, global.objectiveValueHistory); +// SimulationResultsReporter::createPlot("Number of Steps", "Objective value", +// global.objectiveValueHistory); return totalError; } @@ -416,7 +387,7 @@ ReducedModelOptimizer::ReducedModelOptimizer( void ReducedModelOptimizer::initializePatterns( PatternGeometry &fullPattern, PatternGeometry &reducedPattern, - const std::unordered_set &reducedModelExcludedEdges) { + const std::unordered_set &reducedModelExcludedEdges,const int& optimizationParameters) { // fullPattern.setLabel("full_pattern_" + fullPattern.getLabel()); // reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel()); assert(fullPattern.VN() == reducedPattern.VN() && @@ -455,31 +426,64 @@ void ReducedModelOptimizer::initializePatterns( This is not very generic */ computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges); createSimulationMeshes(copyFullPattern, copyReducedPattern); - initializeOptimizationParameters(m_pReducedPatternSimulationMesh); + initializeOptimizationParameters(m_pReducedPatternSimulationMesh,optimizationParameters); } +void updateMesh(long n, const double *x) { + std::shared_ptr &pReducedPatternSimulationMesh = + global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]] + ->pMesh; + + const double E=global.initialParameters(0)*x[0]; + const double A=global.initialParameters(1) * x[1]; + const double beamWidth=std::sqrt(A); + const double beamHeight=beamWidth; + +const double J=global.initialParameters(2) * x[2]; +const double I2=global.initialParameters(3) * x[3]; +const double I3=global.initialParameters(4) * x[4]; + for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { + Element &e = pReducedPatternSimulationMesh->elements[ei]; + e.setDimensions( + RectangularBeamDimensions(beamWidth, + beamHeight)); + e.setMaterial(ElementMaterial(e.material.poissonsRatio, + E)); + e.J = J; + e.I2 = I2; + e.I3 = I3; + } + assert(pReducedPatternSimulationMesh->EN() == 12); + assert(n>=2); + auto R = vcg::RotationMatrix( + ReducedModelOptimizer::patternPlaneNormal, + vcg::math::ToRad(x[n-1] - global.innerHexagonInitialRotationAngle)); + for (int rotationCounter = 0; + rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) { + pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() = + R * global.g_innerHexagonVectors[rotationCounter] * x[n-2]; + } + + pReducedPatternSimulationMesh->reset(); +#ifdef POLYSCOPE_DEFINED + pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); +#endif +} + + void ReducedModelOptimizer::initializeOptimizationParameters( - const std::shared_ptr &mesh) { - global.numberOfOptimizationParameters = 7; + const std::shared_ptr &mesh,const int& optimizationParamters) { + global.numberOfOptimizationParameters = optimizationParamters; global.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; - // } + global.initialParameters(0) = mesh->elements[0].material.youngsModulus; global.initialParameters(1) = mesh->elements[0].A; global.initialParameters(2) = mesh->elements[0].J; global.initialParameters(3) = mesh->elements[0].I2; global.initialParameters(4) = mesh->elements[0].I3; - global.initialParameters(5) = global.innerHexagonInitialPos; + global.initialParameters(optimizationParamters-2) = global.innerHexagonInitialPos; global.innerHexagonInitialRotationAngle = 30; - global.initialParameters(6) = global.innerHexagonInitialRotationAngle; + global.initialParameters(optimizationParamters-1) = global.innerHexagonInitialRotationAngle; } void ReducedModelOptimizer::computeReducedModelSimulationJob( @@ -522,7 +526,7 @@ void ReducedModelOptimizer::visualizeResults( &fullPatternSimulationJobs, const std::vector> &reducedPatternSimulationJobs, - const std::vector &simulationScenarios, + const std::vector &simulationScenarios, const std::unordered_map &reducedToFullInterfaceViMap) { FormFinder simulator; @@ -549,7 +553,7 @@ void ReducedModelOptimizer::visualizeResults( simulator.executeSimulation(pReducedPatternSimulationJob); double normalizationFactor = 1; if (global.optimizationSettings.normalizationStrategy != - Settings::NormalizationStrategy::NonNormalized) { + ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) { normalizationFactor = global.objectiveNormalizationValues[simulationScenarioIndex]; } @@ -559,7 +563,7 @@ void ReducedModelOptimizer::visualizeResults( reducedModelResults.displacements, fullModelResults.displacements, reducedToFullInterfaceViMap, normalizationFactor); std::cout << "Error of simulation scenario " - << simulationScenarioStrings[simulationScenarioIndex] << " is " + << ReducedModelOptimization::simulationScenarioStrings[simulationScenarioIndex] << " is " << error << std::endl; totalError += error; reducedModelResults.registerForDrawing(); @@ -570,7 +574,7 @@ void ReducedModelOptimizer::visualizeResults( "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" "Images/" + pFullPatternSimulationMesh->getLabel() + "_" + - simulationScenarioStrings[simulationScenarioIndex]; + ReducedModelOptimization::simulationScenarioStrings[simulationScenarioIndex]; polyscope::show(); polyscope::screenshot(screenshotFilename, false); fullModelResults.unregister(); @@ -636,10 +640,10 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( } } -ReducedModelOptimizer::Results +ReducedModelOptimization::Results ReducedModelOptimizer::runOptimization(const Settings &settings) { - global.gObjectiveValueHistory.clear(); + global.objectiveValueHistory.clear(); dlib::matrix xMin(global.numberOfOptimizationParameters); dlib::matrix xMax(global.numberOfOptimizationParameters); for (int i = 0; i < global.numberOfOptimizationParameters; i++) { @@ -657,36 +661,43 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) { auto end = std::chrono::system_clock::now(); auto elapsed = std::chrono::duration_cast(end - start); - Results results; + ReducedModelOptimization::Results results; results.numberOfSimulationCrashes = global.numOfSimulationCrashes; - results.x = global.minX; + results.optimalXNameValuePairs.reserve(settings.xRanges.size()); + for(int xVariableIndex=0;xVariableIndex optimalX(results.optimalXNameValuePairs.size()); + for(int xVariableIndex=0;xVariableIndex> fixedVertices; std::unordered_map nodalForces; - const double forceMagnitude = 1; + const double forceMagnitude = 10; //// Axial SimulationScenario scenarioName = SimulationScenario::Axial; @@ -722,7 +733,7 @@ ReducedModelOptimizer::createScenarios( nodalForces[viPair.first] = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) * - forceMagnitude * 2; + forceMagnitude * 10; fixedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; } @@ -756,7 +767,7 @@ ReducedModelOptimizer::createScenarios( nodalForces[viPair.first] = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) * - forceMagnitude * 1; + forceMagnitude * 4; fixedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; } @@ -806,7 +817,7 @@ ReducedModelOptimizer::createScenarios( fixedVertices.clear(); nodalForces.clear(); for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { - nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * 1; + nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * 0.5; fixedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; } @@ -814,7 +825,7 @@ ReducedModelOptimizer::createScenarios( SimulationJob(pMesh, simulationScenarioStrings[scenarioName], fixedVertices, nodalForces, {})); - //// Double using moments + //// Dome using moments scenarioName = SimulationScenario::Dome; fixedVertices.clear(); nodalForces.clear(); @@ -903,7 +914,7 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() { } } -ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( +Results ReducedModelOptimizer::optimize( const Settings &optimizationSettings, const std::vector &simulationScenarios) { diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 400db52..4b3e5e7 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -1,13 +1,14 @@ #ifndef REDUCEDMODELOPTIMIZER_HPP #define REDUCEDMODELOPTIMIZER_HPP -#include "beamformfinder.hpp" +#include "drmsimulationmodel.hpp" #include "csvfile.hpp" #include "edgemesh.hpp" -#include "elementalmesh.hpp" +#include "simulationmesh.hpp" #include "linearsimulationmodel.hpp" #include "matplot/matplot.h" #include +#include "reducedmodeloptimizer_structs.hpp" #ifdef POLYSCOPE_DEFINED #include "polyscope/color_management.h" @@ -16,6 +17,7 @@ using FullPatternVertexIndex = VertexIndex; using ReducedPatternVertexIndex = VertexIndex; class ReducedModelOptimizer { + std::shared_ptr m_pReducedPatternSimulationMesh; std::shared_ptr m_pFullPatternSimulationMesh; std::unordered_map @@ -25,318 +27,14 @@ class ReducedModelOptimizer { std::unordered_map nodeToSlot; std::unordered_map> slotToNode; #ifdef POLYSCOPE_DEFINED - struct StaticColors { - glm::vec3 fullInitial; - glm::vec3 fullDeformed; - glm::vec3 reducedInitial; - glm::vec3 reducedDeformed; - StaticColors() { - fullInitial = {0.416666, 0.109804, 0.890196}; - fullDeformed = {0.583333, 0.890196, 0.109804}; - reducedInitial = {0.890196, 0.109804, 0.193138}; - reducedDeformed = {0.109804, 0.890196, 0.806863}; - } - }; - inline static StaticColors colors; #endif // POLYSCOPE_DEFINED public: inline static int fanSize{6}; inline static VectorType patternPlaneNormal{0, 0, 1}; - enum SimulationScenario { - Axial, - Shear, - Bending, - Dome, - Saddle, - NumberOfSimulationScenarios - }; - struct xRange{ - std::string label; - double min; - double max; - std::string toString() const { - return label + "=[" + std::to_string(min) + "," + std::to_string(max) + - "]"; - } - }; - struct Results; - - struct Settings { - enum NormalizationStrategy { - NonNormalized, - Epsilon, - MaxDisplacement, - EqualDisplacements - }; - inline static vector normalizationStrategyStrings{ - "NonNormalized", "Epsilon", "MaxDsiplacement", "EqualDisplacements"}; - std::vector xRanges; - int numberOfFunctionCalls{100}; - double solutionAccuracy{1e-2}; - NormalizationStrategy normalizationStrategy{Epsilon}; - double normalizationParameter{0.003}; - - std::string toString() const { - std::string settingsString; - if (!xRanges.empty()) { - std::string xRangesString; - for (const xRange &range : xRanges) { - xRangesString += range.toString() + " "; - } - settingsString += xRangesString; - } - settingsString += - "FuncCalls=" + std::to_string(numberOfFunctionCalls) + - " Accuracy=" + std::to_string(solutionAccuracy) + - " Norm=" + normalizationStrategyStrings[normalizationStrategy]; - - return settingsString; - } - - void writeHeaderTo(csvFile &os) const { - if (!xRanges.empty()) { - for (const xRange &range : xRanges) { - os << range.label + " max"; - os << range.label + " min"; - } - } - os << "Function Calls"; - os << "Solution Accuracy"; - os << "Normalization strategy"; - // os << std::endl; - } - - void writeSettingsTo(csvFile &os) const { - if (!xRanges.empty()) { - for (const xRange &range : xRanges) { - os << range.max; - os << range.min; - } - } - os << numberOfFunctionCalls; - os << solutionAccuracy; - os << normalizationStrategyStrings[normalizationStrategy] + "_" + - std::to_string(normalizationParameter); - } - }; - struct Results { - double time{-1}; - int numberOfSimulationCrashes{0}; - std::vector x; - double objectiveValue; - double rawObjectiveValue; - std::vector objectiveValuePerSimulationScenario; - std::vector> fullPatternSimulationJobs; - std::vector> reducedPatternSimulationJobs; - - void save(const string &saveToPath) const { - assert(std::filesystem::is_directory(saveToPath)); - - const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); - for (int simulationJobIndex = 0; - simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationJobIndex]; - std::filesystem::path simulationJobFolderPath( - std::filesystem::path(saveToPath) - .append(pFullPatternSimulationJob->label)); - std::filesystem::create_directory(simulationJobFolderPath); - const auto fullPatternDirectoryPath = - std::filesystem::path(simulationJobFolderPath).append("Full"); - std::filesystem::create_directory(fullPatternDirectoryPath); - pFullPatternSimulationJob->save(fullPatternDirectoryPath.string()); - const std::shared_ptr &pReducedPatternSimulationJob = - reducedPatternSimulationJobs[simulationJobIndex]; - const auto reducedPatternDirectoryPath = - std::filesystem::path(simulationJobFolderPath).append("Reduced"); - if (!std::filesystem::exists(reducedPatternDirectoryPath)) { - std::filesystem::create_directory(reducedPatternDirectoryPath); - } - pReducedPatternSimulationJob->save( - reducedPatternDirectoryPath.string()); - } - } - void load(const string &loadFromPath) { - assert(std::filesystem::is_directory(loadFromPath)); - - for (const auto &directoryEntry : - filesystem::directory_iterator(loadFromPath)) { - const auto simulationScenarioPath = directoryEntry.path(); - if (!std::filesystem::is_directory(simulationScenarioPath)) { - continue; - } - // Load reduced pattern files - for (const auto &fileEntry : filesystem::directory_iterator( - std::filesystem::path(simulationScenarioPath) - .append("Full"))) { - const auto filepath = fileEntry.path(); - if (filepath.extension() == ".json") { - SimulationJob job; - job.load(filepath.string()); - fullPatternSimulationJobs.push_back( - std::make_shared(job)); - } - } - - // Load full pattern files - for (const auto &fileEntry : filesystem::directory_iterator( - std::filesystem::path(simulationScenarioPath) - .append("Reduced"))) { - const auto filepath = fileEntry.path(); - if (filepath.extension() == ".json") { - SimulationJob job; - job.load(filepath.string()); - reducedPatternSimulationJobs.push_back( - std::make_shared(job)); - } - } - } - } -#if POLYSCOPE_DEFINED - 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, false); - - const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); - for (int simulationJobIndex = 0; - simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { - // Drawing of full pattern results - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationJobIndex]; - pFullPatternSimulationJob->registerForDrawing( - fullPatternSimulationJobs[0]->pMesh->getLabel()); - 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 = - "/home/iason/Coding/Projects/Approximating shapes with flat " - "patterns/RodModelOptimizationForPatterns/Results/Images/" + - fullPatternSimulationJobs[0]->pMesh->getLabel() + "_" + - pFullPatternSimulationJob->getLabel(); - 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 " - // << - // simula simulationScenarioStrings[simulationScenarioIndex] - // << " is " - // << error << std::endl; - } - } -#endif // POLYSCOPE_DEFINED - void saveMeshFiles() const { - const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); - assert(numberOfSimulationJobs != 0 && - fullPatternSimulationJobs.size() == - reducedPatternSimulationJobs.size()); - - fullPatternSimulationJobs[0]->pMesh->savePly( - "FullPattern_undeformed.ply"); - reducedPatternSimulationJobs[0]->pMesh->savePly( - "ReducedPattern_undeformed.ply"); - FormFinder simulator; - for (int simulationJobIndex = 0; - simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { - // Drawing of full pattern results - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationJobIndex]; - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob); - fullModelResults.saveDeformedModel(); - - // Drawing of reduced pattern results - const std::shared_ptr &pReducedPatternSimulationJob = - reducedPatternSimulationJobs[simulationJobIndex]; - SimulationResults reducedModelResults = - simulator.executeSimulation(pReducedPatternSimulationJob); - reducedModelResults.saveDeformedModel(); - } - } - - void - writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization, - csvFile &os) { - os << "Total raw Obj value"; - os << "Total Obj value"; - for (int simulationScenarioIndex = 0; - simulationScenarioIndex < - SimulationScenario::NumberOfSimulationScenarios; - simulationScenarioIndex++) { - os << "Obj value " + simulationScenarioStrings[simulationScenarioIndex]; - } - for (const ReducedModelOptimizer::xRange &range : - settings_optimization.xRanges) { - os << range.label; - } - os << "Time(s)"; - os << "#Crashes"; - } - - void - writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization, - csvFile &os) const { - os << rawObjectiveValue; - os << objectiveValue; - for (double scenarioObjValue : objectiveValuePerSimulationScenario) { - os << scenarioObjValue; - } - for (const double &optimalX : x) { - os << optimalX; - } - - for (int unusedXVarCounter = 0; - unusedXVarCounter < settings_optimization.xRanges.size() - x.size(); - unusedXVarCounter++) { - os << "-"; - } - - os << time; - if (numberOfSimulationCrashes == 0) { - os << "-"; - } else { - os << numberOfSimulationCrashes; - } - } - }; - - inline static const std::string simulationScenarioStrings[] = { - "Axial", "Shear", "Bending", "Dome", "Saddle"}; - Results optimize(const Settings &xRanges, - const std::vector &simulationScenarios = - std::vector()); + ReducedModelOptimization::Results optimize(const ReducedModelOptimization::Settings &xRanges, + const std::vector &simulationScenarios = + std::vector()); double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot); @@ -348,9 +46,8 @@ public: SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); - void initializePatterns( - PatternGeometry &fullPattern, PatternGeometry &reducedPatterm, - const std::unordered_set &reducedModelExcludedEges); + void initializePatterns(PatternGeometry &fullPattern, PatternGeometry &reducedPatterm, + const std::unordered_set &reducedModelExcludedEges, const int &optimizationParameters); static void runSimulation(const std::string &filename, std::vector &x); @@ -383,7 +80,7 @@ public: &fullPatternSimulationJobs, const std::vector> &reducedPatternSimulationJobs, - const std::vector &simulationScenarios, + const std::vector &simulationScenarios, const std::unordered_map &reducedToFullInterfaceViMap); @@ -409,12 +106,13 @@ public: &reducedToFullInterfaceViMap, const double &normalizationFactor); -private: + static double objective(double E, double A, double innerHexagonSize, double innerHexagonRotationAngle); + private: static void computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, const std::unordered_map &displacementsReducedToFullMap, Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); - static Results runOptimization(const Settings &settings); + static ReducedModelOptimization::Results runOptimization(const ReducedModelOptimization::Settings &settings); std::vector> createScenarios(const std::shared_ptr &pMesh); void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern, @@ -422,10 +120,11 @@ private: void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); static void - initializeOptimizationParameters(const std::shared_ptr &mesh); + initializeOptimizationParameters(const std::shared_ptr &mesh, const int &optimizationParamters); static double objective(long n, const double *x); FormFinder simulator; void computeObjectiveValueNormalizationFactors(); }; +void updateMesh(long n, const double *x); #endif // REDUCEDMODELOPTIMIZER_HPP