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

Refactoring

This commit is contained in:
iasonmanolas 2021-12-23 15:51:23 +02:00
parent 56e5e043f6
commit c9fc6ccd08
6 changed files with 417 additions and 313 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
CMakeLists.txt
View File

@ -15,28 +15,19 @@ if(NOT EXISTS ${EXTERNAL_DEPS_DIR})
endif() endif()
##Create directory for the external libraries ##Create directory for the external libraries
file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR}) file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR})
#message(${POLYSCOPE_ALREADY_COMPILED})
if(${USE_POLYSCOPE}) if(${USE_POLYSCOPE})
download_project(PROJ POLYSCOPE set(POLYSCOPE_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/polyscope)
GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git download_project(PROJ POLYSCOPE
GIT_TAG master GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git
PREFIX ${EXTERNAL_DEPS_DIR} GIT_TAG master
${UPDATE_DISCONNECTED_IF_AVAILABLE} PREFIX ${EXTERNAL_DEPS_DIR}
) BINARY_DIR ${POLYSCOPE_BINARY_DIR}
add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR}) ${UPDATE_DISCONNECTED_IF_AVAILABLE}
add_compile_definitions(POLYSCOPE_DEFINED) )
add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR})
add_compile_definitions(POLYSCOPE_DEFINED)
endif() endif()
#dlib
set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin)
file(MAKE_DIRECTORY ${DLIB_BIN_DIR})
download_project(PROJ DLIB
GIT_REPOSITORY https://github.com/davisking/dlib.git
GIT_TAG master
BINARY_DIR ${DLIB_BIN_DIR}
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR})
##vcglib devel branch ##vcglib devel branch
download_project(PROJ vcglib_devel download_project(PROJ vcglib_devel
GIT_REPOSITORY https://github.com/IasonManolas/vcglib.git GIT_REPOSITORY https://github.com/IasonManolas/vcglib.git
@ -45,27 +36,36 @@ download_project(PROJ vcglib_devel
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
add_subdirectory(${vcglib_devel_SOURCE_DIR} ${vcglib_devel_BINARY_DIR}) add_subdirectory(${vcglib_devel_SOURCE_DIR} ${vcglib_devel_BINARY_DIR})
##matplot++ lib ##matplot++ lib
set(MATPLOTPLUSPLUS_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/matplot)
download_project(PROJ MATPLOTPLUSPLUS download_project(PROJ MATPLOTPLUSPLUS
GIT_REPOSITORY https://github.com/alandefreitas/matplotplusplus GIT_REPOSITORY https://github.com/alandefreitas/matplotplusplus
GIT_TAG master GIT_TAG master
BINARY_DIR ${MATPLOTPLUSPLUS_BINARY_DIR}
PREFIX ${EXTERNAL_DEPS_DIR} PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
add_subdirectory(${MATPLOTPLUSPLUS_SOURCE_DIR} ${MATPLOTPLUSPLUS_BINARY_DIR}) add_subdirectory(${MATPLOTPLUSPLUS_SOURCE_DIR} ${MATPLOTPLUSPLUS_BINARY_DIR})
##threed-beam-fea ##threed-beam-fea
set(threed-beam-fea_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/threed-beam-fea)
download_project(PROJ threed-beam-fea download_project(PROJ threed-beam-fea
GIT_REPOSITORY https://github.com/IasonManolas/threed-beam-fea.git GIT_REPOSITORY https://github.com/IasonManolas/threed-beam-fea.git
GIT_TAG master GIT_TAG master
BINARY_DIR ${threed-beam-fea_BINARY_DIR}
PREFIX ${EXTERNAL_DEPS_DIR} PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR}) add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR})
##TBB ##TBB
set(TBB_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/tbb)
download_project(PROJ TBB download_project(PROJ TBB
GIT_REPOSITORY https://github.com/wjakob/tbb.git GIT_REPOSITORY https://github.com/wjakob/tbb.git
GIT_TAG master GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR} PREFIX ${EXTERNAL_DEPS_DIR}
BINARY_DIR ${TBB_BINARY_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF) option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF)
@ -83,13 +83,12 @@ target_include_directories(${PROJECT_NAME}
PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph
PUBLIC ${vcglib_devel_SOURCE_DIR} PUBLIC ${vcglib_devel_SOURCE_DIR}
PUBLIC ${threed-beam-fea_SOURCE_DIR} PUBLIC ${threed-beam-fea_SOURCE_DIR}
# PUBLIC ${MySourcesFiles}
) )
if(${MYSOURCES_STATIC_LINK}) if(${MYSOURCES_STATIC_LINK})
message("Linking statically") message("Linking statically")
target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread gfortran quadmath) target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread gfortran quadmath)
else() else()
target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib ThreedBeamFEA tbb pthread) target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot ThreedBeamFEA tbb pthread)
if(${USE_POLYSCOPE}) if(${USE_POLYSCOPE})
target_link_libraries(${PROJECT_NAME} polyscope) target_link_libraries(${PROJECT_NAME} polyscope)
endif() endif()
@ -98,9 +97,11 @@ target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/l
if(USE_ENSMALLEN) if(USE_ENSMALLEN)
##ENSMALLEN ##ENSMALLEN
set(ENSMALLEN_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/ensmallen)
download_project(PROJ ENSMALLEN download_project(PROJ ENSMALLEN
GIT_REPOSITORY https://github.com/mlpack/ensmallen.git GIT_REPOSITORY https://github.com/mlpack/ensmallen.git
GIT_TAG master GIT_TAG master
BINARY_DIR ${ENSMALLEN_BINARY_DIR}
PREFIX ${EXTERNAL_DEPS_DIR} PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
@ -110,7 +111,6 @@ if(USE_ENSMALLEN)
if(${MYSOURCES_STATIC_LINK}) if(${MYSOURCES_STATIC_LINK})
target_link_libraries(${PROJECT_NAME} PUBLIC "-static" ensmallen ${EXTERNAL_DEPS_DIR}/armadillo_build/libarmadillo.a) target_link_libraries(${PROJECT_NAME} PUBLIC "-static" ensmallen ${EXTERNAL_DEPS_DIR}/armadillo_build/libarmadillo.a)
else() else()
target_link_libraries(${PROJECT_NAME} PUBLIC)
target_link_libraries(${PROJECT_NAME} PUBLIC armadillo ensmallen) target_link_libraries(${PROJECT_NAME} PUBLIC armadillo ensmallen)
endif() endif()
target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include) target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include)

5
drmsimulationmodel.cpp
View File

@ -1833,10 +1833,13 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptr<Simul
results.debug_q_normal[vi] = q_normal; results.debug_q_normal[vi] = q_normal;
results.debug_q_nr[vi] = q_nr_nInit; results.debug_q_nr[vi] = q_nr_nInit;
results.rotationalDisplacementQuaternion[vi] results.rotationalDisplacementQuaternion[vi]
//Eigen::Quaterniond R
= (q_normal = (q_normal
* (q_f1_nInit * q_nr_nInit)); //q_f1_nDeformed * q_nr_nDeformed * q_normal; * (q_f1_nInit * q_nr_nInit)); //q_f1_nDeformed * q_nr_nDeformed * q_normal;
//Update the displacement vector to contain the euler angles //Update the displacement vector to contain the euler angles
const Eigen::Vector3d eulerAngles = results.rotationalDisplacementQuaternion[vi] const Eigen::Vector3d eulerAngles = results
.rotationalDisplacementQuaternion[vi]
// R
.toRotationMatrix() .toRotationMatrix()
.eulerAngles(0, 1, 2); .eulerAngles(0, 1, 2);
results.displacements[vi][3] = eulerAngles[0]; results.displacements[vi][3] = eulerAngles[0];

263
reducedmodeloptimizer_structs.hpp
View File

@ -62,140 +62,159 @@ struct xRange
} }
}; };
struct Settings enum OptimizationParameterIndex { E, A, I2, I3, J, R, Theta, NumberOfOptimizationParameters };
{
enum NormalizationStrategy { NonNormalized, Epsilon };
std::vector<xRange> parameterRanges;
inline static vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
int numberOfFunctionCalls{100000};
double solverAccuracy{1e-2};
NormalizationStrategy normalizationStrategy{Epsilon};
double translationNormalizationParameter{0.003};
double rotationNormalizationParameter{3};
bool splitGeometryMaterialOptimization{false};
struct ObjectiveWeights
{
double translational{1};
double rotational{1};
} objectiveWeights;
struct JsonKeys struct Settings
{ {
inline static std::string filename{"OptimizationSettings.json"}; std::filesystem::path intermediateResultsDirectoryPath;
inline static std::string OptimizationVariables{"OptimizationVariables"}; std::vector<std::vector<OptimizationParameterIndex>> optimizationVariables;
inline static std::string NumberOfFunctionCalls{"NumberOfFunctionCalls"}; enum NormalizationStrategy { NonNormalized, Epsilon };
inline static std::string SolverAccuracy{"SolverAccuracy"}; std::vector<xRange> parameterRanges;
inline static std::string TranslationalObjectiveWeight{"TransObjWeight"}; inline static vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
}; int numberOfFunctionCalls{100000};
double solverAccuracy{1e-2};
NormalizationStrategy normalizationStrategy{Epsilon};
double translationNormalizationParameter{0.003};
double rotationNormalizationParameter{3};
struct ObjectiveWeights
{
double translational{1};
double rotational{1};
} objectiveWeights;
void save(const std::filesystem::path &saveToPath) struct JsonKeys
{ {
assert(std::filesystem::is_directory(saveToPath)); inline static std::string filename{"OptimizationSettings.json"};
inline static std::string OptimizationVariables{"OptimizationVariables"};
inline static std::string NumberOfFunctionCalls{"NumberOfFunctionCalls"};
inline static std::string SolverAccuracy{"SolverAccuracy"};
inline static std::string TranslationalObjectiveWeight{"TransObjWeight"};
inline static std::string OptimizationParameters{"OptimizationParameters"};
};
nlohmann::json json; void save(const std::filesystem::path &saveToPath)
//write x ranges {
for (size_t xRangeIndex = 0; xRangeIndex < parameterRanges.size(); xRangeIndex++) { assert(std::filesystem::is_directory(saveToPath));
const auto &xRange = parameterRanges[xRangeIndex];
json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)]
= xRange.toString();
}
json[JsonKeys::NumberOfFunctionCalls] = numberOfFunctionCalls; nlohmann::json json;
json[JsonKeys::SolverAccuracy] = solverAccuracy; //write x ranges
json[JsonKeys::TranslationalObjectiveWeight] = objectiveWeights.translational; for (size_t xRangeIndex = 0; xRangeIndex < parameterRanges.size(); xRangeIndex++) {
const auto &xRange = parameterRanges[xRangeIndex];
json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)]
= xRange.toString();
}
std::filesystem::path jsonFilePath( json[JsonKeys::NumberOfFunctionCalls] = numberOfFunctionCalls;
std::filesystem::path(saveToPath).append(JsonKeys::filename)); json[JsonKeys::SolverAccuracy] = solverAccuracy;
std::ofstream jsonFile(jsonFilePath.string()); json[JsonKeys::TranslationalObjectiveWeight] = objectiveWeights.translational;
jsonFile << json; json[JsonKeys::OptimizationParameters] = optimizationVariables;
}
bool load(const std::filesystem::path &loadFromPath) std::filesystem::path jsonFilePath(
{ std::filesystem::path(saveToPath).append(JsonKeys::filename));
assert(std::filesystem::is_directory(loadFromPath)); std::ofstream jsonFile(jsonFilePath.string());
//Load optimal X jsonFile << json;
nlohmann::json json; }
std::filesystem::path jsonFilepath(
std::filesystem::path(loadFromPath).append(JsonKeys::filename));
if (!std::filesystem::exists(jsonFilepath)) {
return false;
}
std::ifstream ifs(jsonFilepath);
ifs >> json;
//read x ranges bool load(const std::filesystem::path &loadFromPath)
size_t xRangeIndex = 0; {
while (true) { assert(std::filesystem::is_directory(loadFromPath));
const std::string jsonXRangeKey = JsonKeys::OptimizationVariables + "_" //Load optimal X
+ std::to_string(xRangeIndex); nlohmann::json json;
if (!json.contains(jsonXRangeKey)) { std::filesystem::path jsonFilepath(
break; std::filesystem::path(loadFromPath).append(JsonKeys::filename));
} if (!std::filesystem::exists(jsonFilepath)) {
xRange x; return false;
x.fromString(json.at(jsonXRangeKey)); }
parameterRanges.push_back(x); std::ifstream ifs(jsonFilepath);
xRangeIndex++; ifs >> json;
}
numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls); //read x ranges
solverAccuracy = json.at(JsonKeys::SolverAccuracy); size_t xRangeIndex = 0;
objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight); while (true) {
objectiveWeights.rotational = 2 - objectiveWeights.translational; const std::string jsonXRangeKey = JsonKeys::OptimizationVariables + "_"
return true; + std::to_string(xRangeIndex);
} if (!json.contains(jsonXRangeKey)) {
break;
}
xRange x;
x.fromString(json.at(jsonXRangeKey));
parameterRanges.push_back(x);
xRangeIndex++;
}
std::string toString() const optimizationVariables
{ = std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>(
std::string settingsString; (json.at(JsonKeys::OptimizationParameters)));
if (!parameterRanges.empty()) { numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls);
std::string xRangesString; solverAccuracy = json.at(JsonKeys::SolverAccuracy);
for (const xRange &range : parameterRanges) { objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight);
xRangesString += range.toString() + " "; objectiveWeights.rotational = 2 - objectiveWeights.translational;
} return true;
settingsString += xRangesString; }
}
settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls)
+ " Accuracy=" + std::to_string(solverAccuracy)
+ " TransWeight=" + std::to_string(objectiveWeights.translational)
+ " RotWeight=" + std::to_string(objectiveWeights.rotational);
return settingsString; std::string toString() const
} {
std::string settingsString;
if (!parameterRanges.empty()) {
std::string xRangesString;
for (const xRange &range : parameterRanges) {
xRangesString += range.toString() + " ";
}
settingsString += xRangesString;
}
settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls)
+ " Accuracy=" + std::to_string(solverAccuracy)
+ " TransWeight=" + std::to_string(objectiveWeights.translational)
+ " RotWeight=" + std::to_string(objectiveWeights.rotational);
void writeHeaderTo(csvFile &os) const return settingsString;
{ }
if (!parameterRanges.empty()) {
for (const xRange &range : parameterRanges) {
os << range.label + " max";
os << range.label + " min";
}
}
os << "Function Calls";
os << "Solution Accuracy";
os << "Normalization strategy";
os << "Trans weight";
os << "Rot weight";
os << "Splitted geo from mat opt";
// os << std::endl;
}
void writeSettingsTo(csvFile &os) const void writeHeaderTo(csvFile &os) const
{ {
if (!parameterRanges.empty()) { if (!parameterRanges.empty()) {
for (const xRange &range : parameterRanges) { for (const xRange &range : parameterRanges) {
os << range.max; os << range.label + " max";
os << range.min; os << range.label + " min";
} }
} }
os << numberOfFunctionCalls; os << "Function Calls";
os << solverAccuracy; os << "Solution Accuracy";
os << normalizationStrategyStrings[normalizationStrategy] + "_" os << "Normalization strategy";
+ std::to_string(translationNormalizationParameter); os << "Trans weight";
os << objectiveWeights.translational; os << "Rot weight";
os << objectiveWeights.rotational; os << "Optimization parameters";
os << (splitGeometryMaterialOptimization == true ? "true" : "false"); // os << std::endl;
} }
};
void writeSettingsTo(csvFile &os) const
{
if (!parameterRanges.empty()) {
for (const xRange &range : parameterRanges) {
os << range.max;
os << range.min;
}
}
os << numberOfFunctionCalls;
os << solverAccuracy;
os << normalizationStrategyStrings[normalizationStrategy] + "_"
+ std::to_string(translationNormalizationParameter);
os << objectiveWeights.translational;
os << objectiveWeights.rotational;
//export optimization parameters
std::vector<std::vector<int>> vv;
for (const std::vector<OptimizationParameterIndex> &v : optimizationVariables) {
std::vector<int> vi;
vi.reserve(v.size());
for (const OptimizationParameterIndex &parameter : v) {
vi.emplace_back(parameter);
}
vv.push_back(vi);
}
os << Utilities::toString(vv);
}
};
inline bool operator==(const Settings &settings1, const Settings &settings2) noexcept inline bool operator==(const Settings &settings1, const Settings &settings2) noexcept
{ {

319
simulation_structs.hpp
View File

@ -539,7 +539,7 @@ struct SimulationResults
const std::filesystem::path outputFolderPath = outputFolder.empty() const std::filesystem::path outputFolderPath = outputFolder.empty()
? std::filesystem::current_path() ? std::filesystem::current_path()
: std::filesystem::path(outputFolder); : std::filesystem::path(outputFolder);
std::cout << "Saving results to:" << outputFolderPath << std::endl; // std::cout << "Saving results to:" << outputFolderPath << std::endl;
std::filesystem::path simulationJobOutputFolderPath std::filesystem::path simulationJobOutputFolderPath
= std::filesystem::path(outputFolderPath).append("SimulationJob"); = std::filesystem::path(outputFolderPath).append("SimulationJob");
std::filesystem::create_directories(simulationJobOutputFolderPath); std::filesystem::create_directories(simulationJobOutputFolderPath);
@ -569,8 +569,173 @@ struct SimulationResults
void load(const std::filesystem::path &loadFromPath, const std::filesystem::path &loadJobFrom) void load(const std::filesystem::path &loadFromPath, const std::filesystem::path &loadJobFrom)
{ {
//load job
pJob->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json").string()); pJob->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json").string());
load(loadFromPath);
}
void load(const std::filesystem::path &loadFromPath, const std::shared_ptr<SimulationJob> &pJob)
{
this->pJob = pJob;
load(loadFromPath);
}
#ifdef POLYSCOPE_DEFINED
void unregister() const
{
if (!polyscope::hasCurveNetwork(getLabel())) {
std::cerr << "No curve network registered with a name: " << getLabel() << std::endl;
std::cerr << "Nothing to remove." << std::endl;
return;
}
pJob->unregister(getLabel());
polyscope::removeCurveNetwork(getLabel());
}
polyscope::CurveNetwork *registerForDrawing(
const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
const bool &shouldEnable = true,
const double &desiredRadius = 0.001,
// const double &desiredRadius = 0.0001,
const bool &shouldShowFrames = false) const
{
PolyscopeInterface::init();
const std::shared_ptr<SimulationMesh> &mesh = pJob->pMesh;
polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh;
if (!polyscope::hasStructure(getLabel())) {
const auto verts = mesh->getEigenVertices();
const auto edges = mesh->getEigenEdges();
polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(),
verts,
edges);
} else {
polyscopeHandle_deformedEdmeMesh = polyscope::getCurveNetwork(getLabel());
}
polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable);
polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, true);
if (desiredColor.has_value()) {
const glm::vec3 desiredColor_glm(desiredColor.value()[0],
desiredColor.value()[1],
desiredColor.value()[2]);
polyscopeHandle_deformedEdmeMesh->setColor(desiredColor_glm);
}
Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3);
Eigen::MatrixX3d framesX(mesh->VN(), 3);
Eigen::MatrixX3d framesY(mesh->VN(), 3);
Eigen::MatrixX3d framesZ(mesh->VN(), 3);
Eigen::MatrixX3d framesX_initial(mesh->VN(), 3);
Eigen::MatrixX3d framesY_initial(mesh->VN(), 3);
Eigen::MatrixX3d framesZ_initial(mesh->VN(), 3);
// std::unordered_set<int> interfaceNodes{1, 3, 5, 7, 9, 11};
// std::unordered_set<int> interfaceNodes{3, 7, 11, 15, 19, 23};
// std::unordered_set<int> interfaceNodes{};
for (VertexIndex vi = 0; vi < mesh->VN(); vi++) {
const Vector6d &nodalDisplacement = displacements[vi];
nodalDisplacements.row(vi) = Eigen::Vector3d(nodalDisplacement[0],
nodalDisplacement[1],
nodalDisplacement[2]);
// Eigen::Quaternion<double> Rx(Eigen::AngleAxis(nodalDisplacement[2],Eigen::Vector3d(1, 0, 0)));
// Eigen::Quaternion<double> Ry(Eigen::AngleAxis(nodalDisplacement[4],Eigen::Vector3d(0, 1, 0)));
// Eigen::Quaternion<double> Rz(Eigen::AngleAxis(nodalDisplacement[5],Eigen::Vector3d(0, 0, 1)));
// Eigen::Quaternion<double> R=Rx*Ry*Rz;
// if (interfaceNodes.contains(vi)) {
auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1);
auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0);
auto deformedFrameX = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(1, 0, 0);
framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0],
deformedFrameX[1],
deformedFrameX[2]);
framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0],
deformedFrameY[1],
deformedFrameY[2]);
framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0],
deformedNormal[1],
deformedNormal[2]);
framesX_initial.row(vi) = Eigen::Vector3d(1, 0, 0);
framesY_initial.row(vi) = Eigen::Vector3d(0, 1, 0);
framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 1);
// } else {
// framesX.row(vi) = Eigen::Vector3d(0, 0, 0);
// framesY.row(vi) = Eigen::Vector3d(0, 0, 0);
// framesZ.row(vi) = Eigen::Vector3d(0, 0, 0);
// framesX_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
// framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
// framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
// }
}
polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices()
+ nodalDisplacements);
const double frameRadius_default = 0.035;
const double frameLength_default = 0.035;
const bool shouldEnable_default = true;
//if(showFramesOn.contains(vi)){
auto polyscopeHandle_frameX = polyscopeHandle_deformedEdmeMesh
->addNodeVectorQuantity("FrameX", framesX);
polyscopeHandle_frameX->setVectorLengthScale(frameLength_default);
polyscopeHandle_frameX->setVectorRadius(frameRadius_default);
polyscopeHandle_frameX->setVectorColor(
/*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
auto polyscopeHandle_frameY = polyscopeHandle_deformedEdmeMesh
->addNodeVectorQuantity("FrameY", framesY);
polyscopeHandle_frameY->setVectorLengthScale(frameLength_default);
polyscopeHandle_frameY->setVectorRadius(frameRadius_default);
polyscopeHandle_frameY->setVectorColor(
/*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
auto polyscopeHandle_frameZ = polyscopeHandle_deformedEdmeMesh
->addNodeVectorQuantity("FrameZ", framesZ);
polyscopeHandle_frameZ->setVectorLengthScale(frameLength_default);
polyscopeHandle_frameZ->setVectorRadius(frameRadius_default);
polyscopeHandle_frameZ->setVectorColor(
/*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel());
if (!polyscopeHandle_initialMesh) {
polyscopeHandle_initialMesh = mesh->registerForDrawing();
}
// auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh
// ->addNodeVectorQuantity("FrameX", framesX_initial);
// polyscopeHandle_frameX_initial->setVectorLengthScale(frameLength_default);
// polyscopeHandle_frameX_initial->setVectorRadius(frameRadius_default);
// polyscopeHandle_frameX_initial->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
// auto polyscopeHandle_frameY_initial = polyscopeHandle_initialMesh
// ->addNodeVectorQuantity("FrameY", framesY_initial);
// polyscopeHandle_frameY_initial->setVectorLengthScale(frameLength_default);
// polyscopeHandle_frameY_initial->setVectorRadius(frameRadius_default);
// polyscopeHandle_frameY_initial->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
// auto polyscopeHandle_frameZ_initial = polyscopeHandle_initialMesh
// ->addNodeVectorQuantity("FrameZ", framesZ_initial);
// polyscopeHandle_frameZ_initial->setVectorLengthScale(frameLength_default);
// polyscopeHandle_frameZ_initial->setVectorRadius(frameRadius_default);
// polyscopeHandle_frameZ_initial->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
// //}
pJob->registerForDrawing(getLabel());
static bool wasExecuted = false;
if (!wasExecuted) {
std::function<void()> callback = [&]() {
static bool showFrames = shouldShowFrames;
if (ImGui::Checkbox("Show Frames", &showFrames) && showFrames) {
polyscopeHandle_frameX->setEnabled(showFrames);
polyscopeHandle_frameY->setEnabled(showFrames);
polyscopeHandle_frameZ->setEnabled(showFrames);
}
};
PolyscopeInterface::addUserCallback(callback);
wasExecuted = true;
}
return polyscopeHandle_deformedEdmeMesh;
}
#endif
private:
void load(const std::filesystem::path &loadFromPath)
{
converged = true; //assuming it has converged
assert(pJob != nullptr);
//load job
//Use the first .eigenBin file for loading the displacements //Use the first .eigenBin file for loading the displacements
for (auto const &entry : std::filesystem::recursive_directory_iterator(loadFromPath)) { for (auto const &entry : std::filesystem::recursive_directory_iterator(loadFromPath)) {
if (filesystem::is_regular_file(entry) && entry.path().extension() == ".eigenBin") { if (filesystem::is_regular_file(entry) && entry.path().extension() == ".eigenBin") {
@ -591,156 +756,6 @@ struct SimulationResults
Eigen::Vector3d::UnitZ()); Eigen::Vector3d::UnitZ());
} }
} }
#ifdef POLYSCOPE_DEFINED
void unregister() const {
if (!polyscope::hasCurveNetwork(getLabel())) {
std::cerr << "No curve network registered with a name: " << getLabel()
<< std::endl;
std::cerr << "Nothing to remove." << std::endl;
return;
}
pJob->unregister(getLabel());
polyscope::removeCurveNetwork(getLabel());
}
polyscope::CurveNetwork *registerForDrawing(
const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
const bool &shouldEnable = true,
const double &desiredRadius = 0.001,
// const double &desiredRadius = 0.0001,
const bool &shouldShowFrames = false) const
{
PolyscopeInterface::init();
const std::shared_ptr<SimulationMesh> &mesh = pJob->pMesh;
polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh;
if (!polyscope::hasStructure(getLabel())) {
const auto verts = mesh->getEigenVertices();
const auto edges = mesh->getEigenEdges();
polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(),
verts,
edges);
} else {
polyscopeHandle_deformedEdmeMesh = polyscope::getCurveNetwork(getLabel());
}
polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable);
polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, true);
if (desiredColor.has_value()) {
const glm::vec3 desiredColor_glm(desiredColor.value()[0],
desiredColor.value()[1],
desiredColor.value()[2]);
polyscopeHandle_deformedEdmeMesh->setColor(desiredColor_glm);
}
Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3);
Eigen::MatrixX3d framesX(mesh->VN(), 3);
Eigen::MatrixX3d framesY(mesh->VN(), 3);
Eigen::MatrixX3d framesZ(mesh->VN(), 3);
Eigen::MatrixX3d framesX_initial(mesh->VN(), 3);
Eigen::MatrixX3d framesY_initial(mesh->VN(), 3);
Eigen::MatrixX3d framesZ_initial(mesh->VN(), 3);
// std::unordered_set<int> interfaceNodes{1, 3, 5, 7, 9, 11};
// std::unordered_set<int> interfaceNodes{3, 7, 11, 15, 19, 23};
std::unordered_set<int> interfaceNodes{};
for (VertexIndex vi = 0; vi < mesh->VN(); vi++) {
const Vector6d &nodalDisplacement = displacements[vi];
nodalDisplacements.row(vi) = Eigen::Vector3d(nodalDisplacement[0],
nodalDisplacement[1],
nodalDisplacement[2]);
if (interfaceNodes.contains(vi)) {
auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1);
auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0);
auto deformedFrameX = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(1, 0, 0);
framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0],
deformedFrameX[1],
deformedFrameX[2]);
framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0],
deformedFrameY[1],
deformedFrameY[2]);
framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0],
deformedNormal[1],
deformedNormal[2]);
framesX_initial.row(vi) = Eigen::Vector3d(1, 0, 0);
framesY_initial.row(vi) = Eigen::Vector3d(0, 1, 0);
framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 1);
} else {
framesX.row(vi) = Eigen::Vector3d(0, 0, 0);
framesY.row(vi) = Eigen::Vector3d(0, 0, 0);
framesZ.row(vi) = Eigen::Vector3d(0, 0, 0);
framesX_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
}
}
polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices()
+ nodalDisplacements);
const double frameRadius_default = 0.035;
const double frameLength_default = 0.035;
const bool shouldEnable_default = true;
//if(showFramesOn.contains(vi)){
auto polyscopeHandle_frameX = polyscopeHandle_deformedEdmeMesh
->addNodeVectorQuantity("FrameX", framesX);
polyscopeHandle_frameX->setVectorLengthScale(frameLength_default);
polyscopeHandle_frameX->setVectorRadius(frameRadius_default);
polyscopeHandle_frameX->setVectorColor(
/*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
auto polyscopeHandle_frameY = polyscopeHandle_deformedEdmeMesh
->addNodeVectorQuantity("FrameY", framesY);
polyscopeHandle_frameY->setVectorLengthScale(frameLength_default);
polyscopeHandle_frameY->setVectorRadius(frameRadius_default);
polyscopeHandle_frameY->setVectorColor(
/*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
auto polyscopeHandle_frameZ = polyscopeHandle_deformedEdmeMesh
->addNodeVectorQuantity("FrameZ", framesZ);
polyscopeHandle_frameZ->setVectorLengthScale(frameLength_default);
polyscopeHandle_frameZ->setVectorRadius(frameRadius_default);
polyscopeHandle_frameZ->setVectorColor(
/*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel());
if (!polyscopeHandle_initialMesh) {
polyscopeHandle_initialMesh = mesh->registerForDrawing();
}
// auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh
// ->addNodeVectorQuantity("FrameX", framesX_initial);
// polyscopeHandle_frameX_initial->setVectorLengthScale(frameLength_default);
// polyscopeHandle_frameX_initial->setVectorRadius(frameRadius_default);
// polyscopeHandle_frameX_initial->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
// auto polyscopeHandle_frameY_initial = polyscopeHandle_initialMesh
// ->addNodeVectorQuantity("FrameY", framesY_initial);
// polyscopeHandle_frameY_initial->setVectorLengthScale(frameLength_default);
// polyscopeHandle_frameY_initial->setVectorRadius(frameRadius_default);
// polyscopeHandle_frameY_initial->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
// auto polyscopeHandle_frameZ_initial = polyscopeHandle_initialMesh
// ->addNodeVectorQuantity("FrameZ", framesZ_initial);
// polyscopeHandle_frameZ_initial->setVectorLengthScale(frameLength_default);
// polyscopeHandle_frameZ_initial->setVectorRadius(frameRadius_default);
// polyscopeHandle_frameZ_initial->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
// //}
pJob->registerForDrawing(getLabel());
static bool wasExecuted = false;
if (!wasExecuted) {
std::function<void()> callback = [&]() {
static bool showFrames = shouldShowFrames;
if (ImGui::Checkbox("Show Frames", &showFrames) && showFrames) {
polyscopeHandle_frameX->setEnabled(showFrames);
polyscopeHandle_frameY->setEnabled(showFrames);
polyscopeHandle_frameZ->setEnabled(showFrames);
}
};
PolyscopeInterface::addUserCallback(callback);
wasExecuted = true;
}
return polyscopeHandle_deformedEdmeMesh;
}
#endif
}; };
#endif // SIMULATIONHISTORY_HPP #endif // SIMULATIONHISTORY_HPP

10
simulationmesh.cpp
View File

@ -499,13 +499,13 @@ double Element::getMass(const double &materialDensity)
} }
void Element::updateRigidity() { void Element::updateRigidity() {
assert(initialLength != 0); // assert(initialLength != 0);
rigidity.axial = material.youngsModulus * dimensions.A / initialLength; rigidity.axial = material.youngsModulus * dimensions.A / initialLength;
assert(rigidity.axial != 0); // assert(rigidity.axial != 0);
rigidity.torsional = material.G * dimensions.inertia.J / initialLength; rigidity.torsional = material.G * dimensions.inertia.J / initialLength;
assert(rigidity.torsional != 0); // assert(rigidity.torsional != 0);
rigidity.firstBending = 2 * material.youngsModulus * dimensions.inertia.I2 / initialLength; rigidity.firstBending = 2 * material.youngsModulus * dimensions.inertia.I2 / initialLength;
assert(rigidity.firstBending != 0); // assert(rigidity.firstBending != 0);
rigidity.secondBending = 2 * material.youngsModulus * dimensions.inertia.I3 / initialLength; rigidity.secondBending = 2 * material.youngsModulus * dimensions.inertia.I3 / initialLength;
assert(rigidity.secondBending != 0); // assert(rigidity.secondBending != 0);
} }

85
utilities.hpp
View File

@ -140,17 +140,85 @@ struct Vector6d : public std::array<double, 6> {
}; };
namespace Utilities { namespace Utilities {
inline void parseIntegers(const std::string &str, std::vector<size_t> &result) { //inline std::vector<std::string> split(std::string text, char delim)
typedef std::regex_iterator<std::string::const_iterator> re_iterator; //{
typedef re_iterator::value_type re_iterated; // std::string line;
// std::vector<std::string> vec;
// std::stringstream ss(text);
// while (std::getline(ss, line, delim)) {
// vec.push_back(line);
// }
// return vec;
//}
std::regex re("(\\d+)"); inline std::string_view leftTrimSpaces(const std::string_view& str)
{
std::string_view trimmedString=str;
const auto pos(str.find_first_not_of(" \t\n\r\f\v"));
trimmedString.remove_prefix(std::min(pos, trimmedString.length()));
return trimmedString;
}
re_iterator rit(str.begin(), str.end(), re); inline std::string_view rightTrimSpaces(const std::string_view& str)
re_iterator rend; {
std::string_view trimmedString=str;
const auto pos(trimmedString.find_last_not_of(" \t\n\r\f\v"));
trimmedString.remove_suffix(std::min(trimmedString.length() - pos - 1,trimmedString.length()));
return trimmedString;
}
std::transform(rit, rend, std::back_inserter(result), inline std::string_view trimLeftAndRightSpaces(std::string_view str)
[](const re_iterated &it) { return std::stoi(it[1]); }); {
std::string_view trimmedString=str;
trimmedString = leftTrimSpaces(trimmedString);
trimmedString = rightTrimSpaces(trimmedString);
return trimmedString;
}
inline std::vector<std::string> split(const std::string& text, std::string delim)
{
std::vector<std::string> vec;
size_t pos = 0, prevPos = 0;
while (1) {
pos = text.find(delim, prevPos);
if (pos == std::string::npos) {
vec.push_back(text.substr(prevPos));
return vec;
}
vec.push_back(text.substr(prevPos, pos - prevPos));
prevPos = pos + delim.length();
}
}
inline std::string toString(const std::vector<std::vector<int>> &vv)
{
std::string s;
s.append("{");
for (const std::vector<int> &v : vv) {
s.append("{");
for (const int &i : v) {
s.append(std::to_string(i) + ",");
}
s.pop_back();
s.append("}");
}
s.append("}");
return s;
}
inline void parseIntegers(const std::string &str, std::vector<size_t> &result)
{
typedef std::regex_iterator<std::string::const_iterator> re_iterator;
typedef re_iterator::value_type re_iterated;
std::regex re("(\\d+)");
re_iterator rit(str.begin(), str.end(), re);
re_iterator rend;
std::transform(rit, rend, std::back_inserter(result), [](const re_iterated &it) {
return std::stoi(it[1]);
});
} }
inline Eigen::MatrixXd toEigenMatrix(const std::vector<Vector6d> &v) { inline Eigen::MatrixXd toEigenMatrix(const std::vector<Vector6d> &v) {
@ -179,7 +247,6 @@ inline std::vector<Vector6d> fromEigenMatrix(const Eigen::MatrixXd &m)
return v; return v;
} }
// std::string convertToLowercase(const std::string &s) { // std::string convertToLowercase(const std::string &s) {
// std::string lowercase; // std::string lowercase;
// std::transform(s.begin(), s.end(), lowercase.begin(), // std::transform(s.begin(), s.end(), lowercase.begin(),