Moved linear simulation model to MySOurces. Refactored CMakeLists
This commit is contained in:
iasonmanolas
parent
73dda7e33f
commit
015dd6ec45
|
|
@ -17,7 +17,11 @@ set(EXTERNAL_DEPS_DIR "/home/iason/Coding/build/external dependencies")
|
||||||
file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR})
|
file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR})
|
||||||
|
|
||||||
##Polyscope
|
##Polyscope
|
||||||
|
if(${CMAKE_BUILD_TYPE} STREQUAL "Release")
|
||||||
|
set(USE_POLYSCOPE FALSE)
|
||||||
|
else()
|
||||||
set(USE_POLYSCOPE TRUE)
|
set(USE_POLYSCOPE TRUE)
|
||||||
|
endif()
|
||||||
if(${USE_POLYSCOPE})
|
if(${USE_POLYSCOPE})
|
||||||
download_project(PROJ POLYSCOPE
|
download_project(PROJ POLYSCOPE
|
||||||
GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git
|
GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git
|
||||||
|
|
|
||||||
|
|
@ -1,215 +0,0 @@
|
||||||
#ifndef LINEARSIMULATIONMODEL_HPP
|
|
||||||
#define LINEARSIMULATIONMODEL_HPP
|
|
||||||
|
|
||||||
//#include "beam.hpp"
|
|
||||||
#include "simulationresult.hpp"
|
|
||||||
#include "threed_beam_fea.h"
|
|
||||||
#include <filesystem>
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
// struct BeamSimulationProperties {
|
|
||||||
// float crossArea;
|
|
||||||
// float I2;
|
|
||||||
// float I3;
|
|
||||||
// float polarInertia;
|
|
||||||
// float G;
|
|
||||||
// // Properties used by fea
|
|
||||||
// float EA;
|
|
||||||
// float EIz;
|
|
||||||
// float EIy;
|
|
||||||
// float GJ;
|
|
||||||
|
|
||||||
// BeamSimulationProperties(const BeamDimensions &dimensions,
|
|
||||||
// const BeamMaterial &material);
|
|
||||||
//};
|
|
||||||
|
|
||||||
// struct NodalForce {
|
|
||||||
// int index;
|
|
||||||
// int dof;
|
|
||||||
// double magnitude;
|
|
||||||
//};
|
|
||||||
|
|
||||||
// struct SimulationJob {
|
|
||||||
// Eigen::MatrixX3d nodes;
|
|
||||||
// Eigen::MatrixX2i elements;
|
|
||||||
// Eigen::MatrixX3d elementalNormals;
|
|
||||||
// Eigen::VectorXi fixedNodes;
|
|
||||||
// std::vector<NodalForce> nodalForces;
|
|
||||||
// std::vector<BeamDimensions> beamDimensions;
|
|
||||||
// std::vector<BeamMaterial> beamMaterial;
|
|
||||||
//};
|
|
||||||
|
|
||||||
// struct SimulationResults {
|
|
||||||
// std::vector<Eigen::VectorXd> edgeForces; ///< Force values per force
|
|
||||||
// component
|
|
||||||
// ///< #force components x #edges
|
|
||||||
// Eigen::MatrixXd
|
|
||||||
// nodalDisplacements; ///< The displacement of each node #nodes x 3
|
|
||||||
// SimulationResults(const fea::Summary &feaSummary);
|
|
||||||
// SimulationResults() {}
|
|
||||||
//};
|
|
||||||
|
|
||||||
class LinearSimulationModel {
|
|
||||||
public:
|
|
||||||
LinearSimulationModel(){
|
|
||||||
|
|
||||||
}
|
|
||||||
static std::vector<fea::Elem>
|
|
||||||
getFeaElements(const std::shared_ptr<SimulationJob> &job) {
|
|
||||||
const int numberOfEdges = job->pMesh->EN();
|
|
||||||
std::vector<fea::Elem> 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<double> 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<fea::Node>
|
|
||||||
getFeaNodes(const std::shared_ptr<SimulationJob> &job) {
|
|
||||||
const int numberOfNodes = job->pMesh->VN();
|
|
||||||
std::vector<fea::Node> 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<fea::BC>
|
|
||||||
getFeaFixedNodes(const std::shared_ptr<SimulationJob> &job) {
|
|
||||||
std::vector<fea::BC> 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<fea::DOF>(dofIndex), 0));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
return boundaryConditions;
|
|
||||||
}
|
|
||||||
static std::vector<fea::Force>
|
|
||||||
getFeaNodalForces(const std::shared_ptr<SimulationJob> &job) {
|
|
||||||
std::vector<fea::Force> 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<Eigen::VectorXd>(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> &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<fea::Tie> ties;
|
|
||||||
|
|
||||||
// also create an empty list of equations
|
|
||||||
std::vector<fea::Equation> 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
|
|
||||||
|
|
@ -174,8 +174,8 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
|
||||||
|
|
||||||
// run simulations
|
// run simulations
|
||||||
double totalError = 0;
|
double totalError = 0;
|
||||||
// LinearSimulationModel simulator;
|
LinearSimulationModel simulator;
|
||||||
FormFinder simulator;
|
// FormFinder simulator;
|
||||||
for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
|
for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
|
||||||
SimulationResults reducedModelResults = simulator.executeSimulation(
|
SimulationResults reducedModelResults = simulator.executeSimulation(
|
||||||
global.reducedPatternSimulationJobs[simulationScenarioIndex]);
|
global.reducedPatternSimulationJobs[simulationScenarioIndex]);
|
||||||
|
|
@ -765,7 +765,7 @@ ReducedModelOptimizer::createScenarios(
|
||||||
nodalForces[viPair.first] =
|
nodalForces[viPair.first] =
|
||||||
Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0,
|
Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0,
|
||||||
0, 0}) *
|
0, 0}) *
|
||||||
forceMagnitude * 8;
|
forceMagnitude * 20;
|
||||||
fixedVertices[viPair.second] =
|
fixedVertices[viPair.second] =
|
||||||
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
|
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
|
||||||
}
|
}
|
||||||
|
|
@ -799,7 +799,7 @@ ReducedModelOptimizer::createScenarios(
|
||||||
nodalForces[viPair.first] =
|
nodalForces[viPair.first] =
|
||||||
Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0,
|
Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0,
|
||||||
0, 0}) *
|
0, 0}) *
|
||||||
forceMagnitude * 8;
|
forceMagnitude * 5;
|
||||||
fixedVertices[viPair.second] =
|
fixedVertices[viPair.second] =
|
||||||
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
|
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
|
||||||
}
|
}
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue