diff --git a/CMakeLists.txt b/CMakeLists.txt index 0cf0cf4..1229e3e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,7 +17,11 @@ set(EXTERNAL_DEPS_DIR "/home/iason/Coding/build/external dependencies") file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR}) ##Polyscope -set(USE_POLYSCOPE TRUE) +if(${CMAKE_BUILD_TYPE} STREQUAL "Release") + set(USE_POLYSCOPE FALSE) +else() + set(USE_POLYSCOPE TRUE) +endif() if(${USE_POLYSCOPE}) download_project(PROJ POLYSCOPE GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git 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/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index ea6a408..4ce510c 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -174,8 +174,8 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // run simulations double totalError = 0; -// LinearSimulationModel simulator; - FormFinder simulator; + LinearSimulationModel simulator; +// FormFinder simulator; for (const int simulationScenarioIndex : global.simulationScenarioIndices) { SimulationResults reducedModelResults = simulator.executeSimulation( global.reducedPatternSimulationJobs[simulationScenarioIndex]); @@ -765,7 +765,7 @@ ReducedModelOptimizer::createScenarios( nodalForces[viPair.first] = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) * - forceMagnitude * 8; + forceMagnitude * 20; fixedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; } @@ -799,7 +799,7 @@ ReducedModelOptimizer::createScenarios( nodalForces[viPair.first] = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) * - forceMagnitude * 8; + forceMagnitude * 5; fixedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; }