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

Compare commits

base: manolas:master
Branches Tags
manolas:master
manolas:reducedOptWithPatternSimUsingChronos
manolas:fixedBaseScenariosMaxForces
manolas:viscousDamping
..
compare: manolas:viscousDamping
Branches Tags
manolas:reducedOptWithPatternSimUsingChronos
manolas:fixedBaseScenariosMaxForces
manolas:master
manolas:viscousDamping

No commits in common. "master" and "viscousDamping" have entirely different histories.
master ... viscousDam

24 changed files with 1162 additions and 2126 deletions
Show Stats Expand all files Collapse all files

84
CMakeLists.txt
View File

@ -1,10 +1,7 @@
cmake_minimum_required(VERSION 2.8)
project(MySources)
file(GLOB MySourcesFiles ${CMAKE_CURRENT_LIST_DIR}/*.hpp ${CMAKE_CURRENT_LIST_DIR}/*.cpp)
add_library(${PROJECT_NAME} ${MySourcesFiles} )
target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20)
set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
#Download external dependencies NOTE: If the user has one of these libs it shouldn't be downloaded again.
include(${CMAKE_CURRENT_LIST_DIR}/cmake/DownloadProject.cmake)
if (CMAKE_VERSION VERSION_LESS 3.2)
@ -12,114 +9,101 @@ if (CMAKE_VERSION VERSION_LESS 3.2)
else()
set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1")
endif()
set(EXTERNAL_DEPS_DIR "/home/iason/Coding/build/external dependencies")
if(NOT EXISTS ${EXTERNAL_DEPS_DIR})
set(EXTERNAL_DEPS_DIR ${CMAKE_CURRENT_BINARY_DIR})
endif()
##Create directory for the external libraries
file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR})
#message(${POLYSCOPE_ALREADY_COMPILED})
if(${USE_POLYSCOPE})
set(POLYSCOPE_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/polyscope)
download_project(PROJ POLYSCOPE
GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git
GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR}
BINARY_DIR ${POLYSCOPE_BINARY_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
if(NOT EXISTS ${POLYSCOPE_BINARY_DIR})
add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR})
add_compile_definitions(POLYSCOPE_DEFINED)
target_sources(${PROJECT_NAME} PUBLIC ${POLYSCOPE_SOURCE_DIR}/deps/imgui/imgui/misc/cpp/imgui_stdlib.cpp)
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
download_project(PROJ vcglib_devel
GIT_REPOSITORY https://gitea-s2i2s.isti.cnr.it/manolas/vcglib.git
GIT_REPOSITORY https://github.com/IasonManolas/vcglib.git
GIT_TAG devel
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${vcglib_devel_SOURCE_DIR} ${vcglib_devel_BINARY_DIR})
target_sources(${PROJECT_NAME} PUBLIC ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)
##matplot++ lib
set(MATPLOTPLUSPLUS_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/matplot)
download_project(PROJ MATPLOTPLUSPLUS
GIT_REPOSITORY https://github.com/alandefreitas/matplotplusplus
GIT_TAG master
BINARY_DIR ${MATPLOTPLUSPLUS_BINARY_DIR}
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${MATPLOTPLUSPLUS_SOURCE_DIR} ${MATPLOTPLUSPLUS_BINARY_DIR})
##threed-beam-fea
set(threed-beam-fea_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/threed-beam-fea)
download_project(PROJ threed-beam-fea
# GIT_REPOSITORY https://github.com/IasonManolas/threed-beam-fea.git
GIT_REPOSITORY https://gitea-s2i2s.isti.cnr.it/manolas/threed-beam-fea.git
GIT_REPOSITORY https://github.com/IasonManolas/threed-beam-fea.git
GIT_TAG master
BINARY_DIR ${threed-beam-fea_BINARY_DIR}
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR})
##TBB
set(TBB_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/tbb)
download_project(PROJ TBB
GIT_REPOSITORY https://github.com/wjakob/tbb.git
GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR}
BINARY_DIR ${TBB_BINARY_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF)
add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR})
link_directories(${TBB_BINARY_DIR})
###Eigen 3 NOTE: Eigen is required on the system the code is ran
find_package(Eigen3 3.3 REQUIRED)
if(MSVC)
add_compile_definitions(_HAS_STD_BYTE=0)
endif(MSVC)
#link_directories(${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
if(${MYSOURCES_STATIC_LINK})
message("Linking statically")
target_link_libraries(${PROJECT_NAME} PUBLIC -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a #[[tbb_static]] pthread gfortran quadmath)
else()
# set_target_properties(${PROJECT_NAME} PROPERTIES POSITION_INDEPENDENT_CODE TRUE)
target_link_libraries(${PROJECT_NAME} PUBLIC Eigen3::Eigen matplot ThreedBeamFEA tbb pthread)
if(${USE_POLYSCOPE})
message("Using polyscope")
target_link_libraries(${PROJECT_NAME} PUBLIC polyscope)
endif()
endif()
target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
#link_directories(${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
file(GLOB MySourcesFiles ${CMAKE_CURRENT_LIST_DIR}/*.hpp ${CMAKE_CURRENT_LIST_DIR}/*.cpp)
add_library(${PROJECT_NAME} ${MySourcesFiles} ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)
target_include_directories(${PROJECT_NAME}
PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph
PUBLIC ${vcglib_devel_SOURCE_DIR}
PUBLIC ${threed-beam-fea_SOURCE_DIR}
PUBLIC matplot
PUBLIC ${MySourcesFiles}
)
if(USE_ENSMALLEN)
##ENSMALLEN
set(ENSMALLEN_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/ensmallen)
download_project(PROJ ENSMALLEN
GIT_REPOSITORY https://github.com/mlpack/ensmallen.git
GIT_TAG master
BINARY_DIR ${ENSMALLEN_BINARY_DIR}
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${ENSMALLEN_SOURCE_DIR} ${ENSMALLEN_BINARY_DIR})
set(ARMADILLO_SOURCE_DIR "/home/iason/Coding/Libraries/armadillo")
add_subdirectory(${ARMADILLO_SOURCE_DIR} ${EXTERNAL_DEPS_DIR}/armadillo_build)
if(${MYSOURCES_STATIC_LINK})
target_link_libraries(${PROJECT_NAME} PUBLIC "-static" ensmallen ${EXTERNAL_DEPS_DIR}/armadillo_build/libarmadillo.a)
target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread)
elseif(${USE_POLYSCOPE})
target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib glad ThreedBeamFEA tbb pthread polyscope)
else()
target_link_libraries(${PROJECT_NAME} PUBLIC armadillo ensmallen)
endif()
target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include ${ENSMALLEN_SOURCE_DIR})
target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib glad ThreedBeamFEA tbb pthread)
endif()
target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)

32
beam.hpp
View File

@ -9,36 +9,15 @@ struct RectangularBeamDimensions {
inline static std::string name{"Rectangular"};
double b;
double h;
double A{0}; // cross sectional area
struct MomentsOfInertia
{
double I2{0}; // second moment of inertia
double I3{0}; // third moment of inertia
double J{0}; // torsional constant (polar moment of inertia)
} inertia;
RectangularBeamDimensions(const double &width, const double &height) : b(width), h(height)
{
RectangularBeamDimensions(const double &width, const double &height)
: b(width), h(height) {
assert(width > 0 && height > 0);
updateProperties();
}
RectangularBeamDimensions() : b(0.002), h(0.002) { updateProperties(); }
RectangularBeamDimensions() : b(0.002), h(0.002) {}
std::string toString() const {
return std::string("b=") + std::to_string(b) + std::string(" h=") +
std::to_string(h);
}
void updateProperties()
{
A = b * h;
inertia.I2 = b * std::pow(h, 3) / 12;
inertia.I3 = h * std::pow(b, 3) / 12;
inertia.J = inertia.I2 + inertia.I3;
}
static void computeMomentsOfInertia(const RectangularBeamDimensions &dimensions,
MomentsOfInertia &inertia);
};
struct CylindricalBeamDimensions {
@ -60,21 +39,18 @@ struct ElementMaterial
{
double poissonsRatio;
double youngsModulus;
double G;
ElementMaterial(const double &poissonsRatio, const double &youngsModulus)
: poissonsRatio(poissonsRatio), youngsModulus(youngsModulus)
{
assert(poissonsRatio <= 0.5 && poissonsRatio >= -1);
updateProperties();
}
ElementMaterial() : poissonsRatio(0.3), youngsModulus(200 * 1e9) { updateProperties(); }
ElementMaterial() : poissonsRatio(0.3), youngsModulus(200 * 1e9) {}
std::string toString() const
{
return std::string("Material:") + std::string("\nPoisson's ratio=")
+ std::to_string(poissonsRatio) + std::string("\nYoung's Modulus(GPa)=")
+ std::to_string(youngsModulus / 1e9);
}
void updateProperties() { G = youngsModulus / (2 * (1 + poissonsRatio)); }
};
#endif // BEAM_HPP

9
csvfile.hpp
View File

@ -22,18 +22,17 @@ class csvFile {
const std::string special_chars_;
public:
csvFile(const std::filesystem::path &filename,
const bool &overwrite,
csvFile(const std::string filename, const bool &overwrite,
const std::string separator = ",")
: fs_(), is_first_(true), separator_(separator), escape_seq_("\""), special_chars_("\"")
{
: fs_(), is_first_(true), separator_(separator), escape_seq_("\""),
special_chars_("\"") {
fs_.exceptions(std::ios::failbit | std::ios::badbit);
if (filename.empty()) {
fs_.copyfmt(std::cout);
fs_.clear(std::cout.rdstate());
fs_.basic_ios<char>::rdbuf(std::cout.rdbuf());
} else {
if (!std::filesystem::exists(filename)) {
if (!std::filesystem::exists(std::filesystem::path(filename))) {
std::ofstream outfile(filename);
outfile.close();
}

256
drmsimulationmodel.cpp
View File

@ -217,16 +217,28 @@ void DRMSimulationModel::runUnitTests()
// polyscope::show();
}
void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob)
void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings)
{
//#ifdef USE_ENSMALLEN
// this->pJob = pJob;
//#endif
mSettings = settings;
mCurrentSimulationStep = 0;
history.clear();
history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
constrainedVertices.clear();
pMesh.reset();
pMesh = std::make_unique<SimulationMesh>(*pJob->pMesh);
vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
plotYValues.clear();
plotHandle.reset();
checkedForMaximumMoment = false;
externalMomentsNorm = 0;
Dt = settings.Dtini;
numOfDampings = 0;
shouldTemporarilyDampForces = false;
externalLoadStep = 1;
isVertexConstrained.clear();
minTotalResidualForcesNorm = std::numeric_limits<double>::max();
constrainedVertices.clear();
constrainedVertices = pJob->constrainedVertices;
if (!pJob->nodalForcedDisplacements.empty()) {
for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair : pJob->nodalForcedDisplacements) {
@ -239,26 +251,6 @@ void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob)
for (auto fixedVertex : constrainedVertices) {
isVertexConstrained[fixedVertex.first] = true;
}
}
void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings)
{
mSettings = settings;
mCurrentSimulationStep = 0;
history.clear();
history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
plotYValues.clear();
plotHandle.reset();
checkedForMaximumMoment = false;
externalMomentsNorm = 0;
Dt = settings.Dtini;
numOfDampings = 0;
shouldTemporarilyDampForces = false;
externalLoadStep = 1;
isVertexConstrained.clear();
minTotalResidualForcesNorm = std::numeric_limits<double>::max();
reset(pJob);
#ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) {
@ -796,21 +788,21 @@ double DRMSimulationModel::computeTotalInternalPotentialEnergy()
const EdgeIndex ei = pMesh->getIndex(e);
const double e_k = element.length - element.initialLength;
const double axialPE = pow(e_k, 2) * element.material.youngsModulus * element.dimensions.A
const double axialPE = pow(e_k, 2) * element.material.youngsModulus * element.A
/ (2 * element.initialLength);
const double theta1Diff = theta1_jplus1 - theta1_j;
const double torsionalPE = element.material.G * element.dimensions.inertia.J
* pow(theta1Diff, 2) / (2 * element.initialLength);
const double torsionalPE = element.G * element.inertia.J * pow(theta1Diff, 2)
/ (2 * element.initialLength);
const double firstBendingPE_inBracketsTerm = 4 * pow(theta2_j, 2)
+ 4 * theta2_j * theta2_jplus1
+ 4 * pow(theta2_jplus1, 2);
const double firstBendingPE = firstBendingPE_inBracketsTerm * element.material.youngsModulus
* element.dimensions.inertia.I2 / (2 * element.initialLength);
* element.inertia.I2 / (2 * element.initialLength);
const double secondBendingPE_inBracketsTerm = 4 * pow(theta3_j, 2)
+ 4 * theta3_j * theta3_jplus1
+ 4 * pow(theta3_jplus1, 2);
const double secondBendingPE = secondBendingPE_inBracketsTerm
* element.material.youngsModulus * element.dimensions.inertia.I3
* element.material.youngsModulus * element.inertia.I3
/ (2 * element.initialLength);
totalInternalPotentialEnergy += axialPE + torsionalPE + firstBendingPE + secondBendingPE;
@ -901,8 +893,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
// #pragma omp parallel for schedule(static) num_threads(6)
for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi]
&& fixedVertices.at(edgeNode.vi)
.contains(dofi);
&& fixedVertices.at(edgeNode.vi).contains(dofi);
if (!isDofConstrainedFor_ev) {
DifferentiateWithRespectTo dui{ev, dofi};
// Axial force computation
@ -1073,6 +1064,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
}
#endif
}
pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm;
pMesh->totalResidualForcesNorm = totalResidualForcesNorm;
pMesh->averageResidualForcesNorm = totalResidualForcesNorm / pMesh->VN();
// pMesh->averageResidualForcesNorm = sumOfResidualForces.norm() / pMesh->VN();
@ -1261,15 +1253,14 @@ void DRMSimulationModel::updateNodalMasses()
for (const EdgePointer &ep : pMesh->nodes[v].incidentElements) {
const size_t ei = pMesh->getIndex(ep);
const Element &elem = pMesh->elements[ep];
const double SkTranslational = elem.material.youngsModulus * elem.dimensions.A
/ elem.length;
const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length;
translationalSumSk += SkTranslational;
const double lengthToThe3 = std::pow(elem.length, 3);
const long double SkRotational_I2 = elem.material.youngsModulus * elem.dimensions.inertia.I2
const long double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2
/ lengthToThe3;
const long double SkRotational_I3 = elem.material.youngsModulus * elem.dimensions.inertia.I3
const long double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3
/ lengthToThe3;
const long double SkRotational_J = elem.material.youngsModulus * elem.dimensions.inertia.J
const long double SkRotational_J = elem.material.youngsModulus * elem.inertia.J
/ lengthToThe3;
rotationalSumSk_I2 += SkRotational_I2;
rotationalSumSk_I3 += SkRotational_I3;
@ -1628,15 +1619,14 @@ void DRMSimulationModel::updatePositionsOnTheFly(
double rotationalSumSk_J = 0;
for (const EdgePointer &ep : pMesh->nodes[v].incidentElements) {
const Element &elem = pMesh->elements[ep];
const double SkTranslational = elem.material.youngsModulus * elem.dimensions.A
/ elem.length;
const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length;
translationalSumSk += SkTranslational;
const double lengthToThe3 = std::pow(elem.length, 3);
const double SkRotational_I2 = elem.material.youngsModulus * elem.dimensions.inertia.I2
const double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const double SkRotational_I3 = elem.material.youngsModulus * elem.dimensions.inertia.I3
const double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const double SkRotational_J = elem.material.youngsModulus * elem.dimensions.inertia.J
const double SkRotational_J = elem.material.youngsModulus * elem.inertia.J
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
rotationalSumSk_I2 += SkRotational_I2;
rotationalSumSk_I3 += SkRotational_I3;
@ -1833,13 +1823,10 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptr<Simul
results.debug_q_normal[vi] = q_normal;
results.debug_q_nr[vi] = q_nr_nInit;
results.rotationalDisplacementQuaternion[vi]
//Eigen::Quaterniond R
= (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
const Eigen::Vector3d eulerAngles = results
.rotationalDisplacementQuaternion[vi]
// R
const Eigen::Vector3d eulerAngles = results.rotationalDisplacementQuaternion[vi]
.toRotationMatrix()
.eulerAngles(0, 1, 2);
results.displacements[vi][3] = eulerAngles[0];
@ -1889,7 +1876,7 @@ void DRMSimulationModel::printCurrentState() const
void DRMSimulationModel::printDebugInfo() const
{
std::cout << pMesh->elements[0].rigidity.toString() << std::endl;
std::cout << "Number of dampings:" << numOfDampings << std::endl;
std::cout << "Number of dampings:" << numOfDampings << endl;
printCurrentState();
}
@ -2012,10 +1999,10 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
std::vector<double> I3(pMesh->EN());
for (const EdgeType &e : pMesh->edge) {
const size_t ei = pMesh->getIndex(e);
A[ei] = pMesh->elements[e].dimensions.A;
J[ei] = pMesh->elements[e].dimensions.inertia.J;
I2[ei] = pMesh->elements[e].dimensions.inertia.I2;
I3[ei] = pMesh->elements[e].dimensions.inertia.I3;
A[ei] = pMesh->elements[e].A;
J[ei] = pMesh->elements[e].inertia.J;
I2[ei] = pMesh->elements[e].inertia.I2;
I3[ei] = pMesh->elements[e].inertia.I3;
}
polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A);
@ -2669,119 +2656,54 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
return results;
}
//#ifdef USE_ENSMALLEN
//void DRMSimulationModel::setJob(const std::shared_ptr<SimulationJob> &pJob)
//{
// reset(pJob);
// updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
// // arma::mat externalForces(pMesh->VN() * NumDoF);
// // for (int vi = 0; vi < pMesh->VN(); vi++) {
// // const Vector6d &nodeExternalForce = pMesh.vert[vi].node.externalForce;
// // for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
// // externalForces(vi * NumDoF + dofi) = nodeExternalForce[dofi];
// // }
// // }
//}
//SimulationMesh *DRMSimulationModel::getDeformedMesh(const arma::mat &x)
//{
// for (int vi = 0; vi < pMesh->VN(); vi++) {
// Node &node = pMesh->nodes[vi];
// for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
// node.displacements[dofi] = x(vi * DoF::NumDoF + dofi);
// }
// }
// applyDisplacements(constrainedVertices);
// if (!pJob->nodalForcedDisplacements.empty()) {
// applyForcedDisplacements(pJob->nodalForcedDisplacements);
// }
// return pMesh.get();
//}
//double DRMSimulationModel::EvaluateWithGradient(const arma::mat &x, arma::mat &g)
//{
// for (int vi = 0; vi < pMesh->VN(); vi++) {
// Node &node = pMesh->nodes[vi];
// for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
// node.displacements[dofi] = x(vi * DoF::NumDoF + dofi);
// }
// }
// applyDisplacements(constrainedVertices);
// if (!pJob->nodalForcedDisplacements.empty()) {
// applyForcedDisplacements(pJob->nodalForcedDisplacements);
// }
// updateElementalLengths();
// updateNormalDerivatives();
// updateT1Derivatives();
// updateRDerivatives();
// updateT2Derivatives();
// updateT3Derivatives();
// updateResidualForcesOnTheFly(constrainedVertices);
// for (int vi = 0; vi < pMesh->VN(); vi++) {
// Node &node = pMesh->nodes[vi];
// for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
// g(vi * NumDoF + dofi) = node.force.residual[dofi];
// }
// }
// mCurrentSimulationStep++;
// const double PE = computeTotalPotentialEnergy();
// std::cout << "PE:" << PE << std::endl;
// return PE;
//}
//#endif
void DRMSimulationModel::Settings::save(const std::filesystem::path &jsonFilePath) const
void DRMSimulationModel::Settings::save(const std::filesystem::path &folderPath) const
{
std::filesystem::create_directories(folderPath);
nlohmann::json json;
json[GET_VARIABLE_NAME(shouldDraw)] = shouldDraw;
json[GET_VARIABLE_NAME(beVerbose)] = beVerbose;
json[GET_VARIABLE_NAME(shouldCreatePlots)] = shouldCreatePlots;
json[GET_VARIABLE_NAME(Dtini)] = Dtini;
json[GET_VARIABLE_NAME(xi)] = xi;
json[GET_VARIABLE_NAME(gamma)] = gamma;
json[GET_VARIABLE_NAME(useKineticDamping)] = totalResidualForcesNormThreshold;
json[jsonLabels.gamma] = gamma;
json[jsonLabels.shouldDraw] = shouldDraw ? "true" : "false";
json[jsonLabels.shouldCreatePlots] = shouldCreatePlots ? "true" : "false";
json[jsonLabels.beVerbose] = beVerbose ? "true" : "false";
json[jsonLabels.Dtini] = Dtini;
json[jsonLabels.xi] = xi;
if (maxDRMIterations.has_value()) {
json[GET_VARIABLE_NAME(jsonLabels.maxDRMIterations)] = maxDRMIterations.value();
json[jsonLabels.maxDRMIterations] = maxDRMIterations.value();
}
if (debugModeStep.has_value()) {
json[GET_VARIABLE_NAME(debugModeStep)] = debugModeStep.value();
json[jsonLabels.debugModeStep] = debugModeStep.value();
}
if (totalTranslationalKineticEnergyThreshold.has_value()) {
json[GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold)]
json[jsonLabels.totalTranslationalKineticEnergyThreshold]
= totalTranslationalKineticEnergyThreshold.value();
}
if (averageResidualForcesCriterionThreshold.has_value()) {
json[GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold)]
json[jsonLabels.averageResidualForcesCriterionThreshold]
= averageResidualForcesCriterionThreshold.value();
}
if (averageResidualForcesCriterionThreshold.has_value()) {
json[jsonLabels.averageResidualForcesCriterionThreshold]
= averageResidualForcesCriterionThreshold.value();
}
if (linearGuessForceScaleFactor.has_value()) {
json[GET_VARIABLE_NAME(linearGuessForceScaleFactor)] = linearGuessForceScaleFactor.value();
json[jsonLabels.linearGuessForceScaleFactor] = linearGuessForceScaleFactor.value();
}
if (viscousDampingFactor.has_value()) {
json[GET_VARIABLE_NAME(viscousDampingFactor)] = viscousDampingFactor.value();
}
json[GET_VARIABLE_NAME(useKineticDamping)] = useKineticDamping;
// if (intermediateResultsSaveStep.has_value()) {
// json[GET_VARIABLE_NAME(intermediateResultsSaveStep)] = intermediateResultsSaveStep.value();
// }
// if (saveIntermediateBestStates.has_value()) {
// json[GET_VARIABLE_NAME(saveIntermediateBestStates)] = saveIntermediateBestStates.value()
// ? "true"
// : "false";
// }
if (saveIntermediateBestStates.has_value()) {
json[GET_VARIABLE_NAME(saveIntermediateBestStates)] = saveIntermediateBestStates.value()
? "true"
: "false";
}
std::filesystem::create_directories(jsonFilePath.parent_path());
std::ofstream jsonFile(jsonFilePath);
std::cout << "Saving DRM settings to:" << jsonFilePath << std::endl;
const std::string jsonFilename = "drmSettings.json";
std::ofstream jsonFile(std::filesystem::path(folderPath).append(jsonFilename));
jsonFile << json;
jsonFile.close();
}
bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePath)
@ -2805,40 +2727,41 @@ bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePat
std::ifstream ifs(jsonFilePath.string());
ifs >> json;
if (json.contains(GET_VARIABLE_NAME(shouldDraw))) {
shouldDraw = json.at(GET_VARIABLE_NAME(shouldDraw)) == "true" ? true : false;
if (json.contains(jsonLabels.shouldDraw)) {
shouldDraw = json.at(jsonLabels.shouldDraw) == "true" ? true : false;
}
if (json.contains(GET_VARIABLE_NAME(beVerbose))) {
beVerbose = json.at(GET_VARIABLE_NAME(beVerbose)) == "true" ? true : false;
if (json.contains(jsonLabels.beVerbose)) {
beVerbose = json.at(jsonLabels.beVerbose) == "true" ? true : false;
}
if (json.contains(GET_VARIABLE_NAME(shouldCreatePlots))) {
shouldCreatePlots = json.at(GET_VARIABLE_NAME(shouldCreatePlots)) == "true" ? true : false;
if (json.contains(jsonLabels.shouldCreatePlots)) {
shouldCreatePlots = json.at(jsonLabels.shouldCreatePlots) == "true" ? true : false;
}
if (json.contains(GET_VARIABLE_NAME(Dtini))) {
Dtini = json.at(GET_VARIABLE_NAME(Dtini));
if (json.contains(jsonLabels.Dtini)) {
Dtini = json.at(jsonLabels.Dtini);
}
if (json.contains(GET_VARIABLE_NAME(xi))) {
xi = json.at(GET_VARIABLE_NAME(xi));
if (json.contains(jsonLabels.xi)) {
xi = json.at(jsonLabels.xi);
}
if (json.contains(GET_VARIABLE_NAME(maxDRMIterations))) {
maxDRMIterations = json.at(GET_VARIABLE_NAME(maxDRMIterations));
if (json.contains(jsonLabels.maxDRMIterations)) {
maxDRMIterations = json.at(jsonLabels.maxDRMIterations);
}
if (json.contains(GET_VARIABLE_NAME(debugModeStep))) {
debugModeStep = json.at(GET_VARIABLE_NAME(debugModeStep));
if (json.contains(jsonLabels.debugModeStep)) {
debugModeStep = json.at(jsonLabels.debugModeStep);
}
if (json.contains(GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold))) {
if (json.contains(jsonLabels.totalTranslationalKineticEnergyThreshold)) {
totalTranslationalKineticEnergyThreshold = json.at(
GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold));
jsonLabels.totalTranslationalKineticEnergyThreshold);
}
if (json.contains(GET_VARIABLE_NAME(gamma))) {
gamma = json.at(GET_VARIABLE_NAME(gamma));
if (json.contains(jsonLabels.gamma)) {
gamma = json.at(jsonLabels.gamma);
}
if (json.contains(GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold))) {
if (json.contains(jsonLabels.averageResidualForcesCriterionThreshold)) {
averageResidualForcesCriterionThreshold = json.at(
GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold));
jsonLabels.averageResidualForcesCriterionThreshold);
}
if (json.contains(GET_VARIABLE_NAME(linearGuessForceScaleFactor))) {
linearGuessForceScaleFactor = json.at(GET_VARIABLE_NAME(linearGuessForceScaleFactor));
if (json.contains(jsonLabels.linearGuessForceScaleFactor)) {
linearGuessForceScaleFactor = json.at(jsonLabels.linearGuessForceScaleFactor);
}
// if (json.contains(GET_VARIABLE_NAME(intermediateResultsSaveStep))) {
@ -2851,7 +2774,6 @@ bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePat
? true
: false;
}
// std::cout << json.dump() << std::endl;
return true;
}

39
drmsimulationmodel.hpp
View File

@ -1,21 +1,13 @@
#ifndef BEAMFORMFINDER_HPP
#define BEAMFORMFINDER_HPP
//#ifdef USE_MATPLOT
//#include "matplot.h"
#include <matplot/matplot.h>
//#endif
#include "matplot/matplot.h"
#include "simulationmesh.hpp"
#include "simulationmodel.hpp"
#include <Eigen/Dense>
#include <filesystem>
#include <unordered_set>
#ifdef USE_ENSMALLEN
#include <armadillo>
#include <ensmallen.hpp>
#endif
struct SimulationJob;
enum DoF { Ux = 0, Uy, Uz, Nx, Ny, Nr, NumDoF };
@ -33,6 +25,7 @@ class DRMSimulationModel : public SimulationModel
public:
struct Settings
{
// bool isDebugMode{false};
bool shouldDraw{false};
bool beVerbose{false};
bool shouldCreatePlots{false};
@ -45,8 +38,8 @@ public:
int gradualForcedDisplacementSteps{50};
// int desiredGradualExternalLoadsSteps{1};
double gamma{0.8};
double totalResidualForcesNormThreshold{1e-8};
double totalExternalForcesNormPercentageTermination{1e-8};
double totalResidualForcesNormThreshold{1e-3};
double totalExternalForcesNormPercentageTermination{1e-3};
std::optional<int> maxDRMIterations;
std::optional<int> debugModeStep;
std::optional<double> totalTranslationalKineticEnergyThreshold;
@ -54,11 +47,11 @@ public:
std::optional<double> linearGuessForceScaleFactor;
// std::optional<int> intermediateResultsSaveStep;
std::optional<bool> saveIntermediateBestStates;
Settings() {}
void save(const std::filesystem::path &folderPath) const;
bool load(const std::filesystem::path &filePath);
std::optional<double> viscousDampingFactor;
bool useKineticDamping{true};
Settings() {}
void save(const std::filesystem::path &jsonFilePath) const;
bool load(const std::filesystem::path &filePath);
struct JsonLabels
{
const std::string shouldDraw{"shouldDraw"};
@ -66,16 +59,14 @@ public:
const std::string shouldCreatePlots{"shouldCreatePlots"};
const std::string Dtini{"DtIni"};
const std::string xi{"xi"};
const std::string gamma{"gamma"};
const std::string totalResidualForcesNormThreshold;
const std::string maxDRMIterations{"maxIterations"};
const std::string debugModeStep{"debugModeStep"};
const std::string totalTranslationalKineticEnergyThreshold{
"totalTranslationaKineticEnergyThreshold"};
const std::string gamma{"gamma"};
const std::string averageResidualForcesCriterionThreshold{
"averageResidualForcesThreshold"};
const std::string linearGuessForceScaleFactor{"linearGuessForceScaleFactor"};
const std::string viscousDampingFactor{"viscousDampingFactor"};
};
static JsonLabels jsonLabels;
};
@ -103,9 +94,7 @@ private:
// Eigen::Tensor<double, 4> theta3Derivatives;
// std::unordered_map<MyKeyType, double, key_hash> theta3Derivatives;
bool shouldApplyInitialDistortion = false;
//#ifdef USE_ENSMALLEN
// std::shared_ptr<SimulationJob> pJob;
//#endif
void reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings);
void updateNodalInternalForces(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
@ -153,7 +142,7 @@ private:
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
#ifdef POLYSCOPE_DEFINED
void draw(const std::string &screenshotsFolder = {});
void draw(const string &screenshotsFolder= {});
#endif
void
updateNodalInternalForce(Vector6d &nodalInternalForce,
@ -252,19 +241,11 @@ private:
virtual SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob);
void reset(const std::shared_ptr<SimulationJob> &pJob);
public:
DRMSimulationModel();
SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings,
const SimulationResults &solutionGuess = SimulationResults());
#ifdef USE_ENSMALLEN
double EvaluateWithGradient(const arma::mat &x, arma::mat &g);
void setJob(const std::shared_ptr<SimulationJob> &pJob);
SimulationMesh *getDeformedMesh(const arma::mat &x);
#endif
static void runUnitTests();
};

57
edgemesh.cpp
View File

@ -25,13 +25,13 @@ Eigen::MatrixX3d VCGEdgeMesh::getEigenEdgeNormals() const {
return eigenEdgeNormals;
}
bool VCGEdgeMesh::save(const std::filesystem::__cxx11::path &meshFilePath)
bool VCGEdgeMesh::save(const string &plyFilename)
{
std::string filename = meshFilePath;
std::string filename = plyFilename;
if (filename.empty()) {
filename = std::filesystem::current_path().append(getLabel() + ".ply").string();
} else if (std::filesystem::is_directory(std::filesystem::path(meshFilePath))) {
filename = std::filesystem::path(meshFilePath).append(getLabel() + ".ply").string();
} else if (std::filesystem::is_directory(std::filesystem::path(plyFilename))) {
filename = std::filesystem::path(plyFilename).append(getLabel() + ".ply").string();
}
assert(std::filesystem::path(filename).extension().string() == ".ply");
unsigned int mask = 0;
@ -197,14 +197,14 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth,
return true;
}
bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath)
{
std::string usedPath = meshFilePath;
if (std::filesystem::path(meshFilePath).is_relative()) {
usedPath = std::filesystem::absolute(meshFilePath).string();
bool VCGEdgeMesh::load(const string &plyFilename) {
std::string usedPath = plyFilename;
if (std::filesystem::path(plyFilename).is_relative()) {
usedPath = std::filesystem::absolute(plyFilename).string();
}
assert(std::filesystem::exists(usedPath));
Clear();
this->Clear();
// if (!loadUsingNanoply(usedPath)) {
// std::cerr << "Error: Unable to open " + usedPath << std::endl;
// return false;
@ -217,15 +217,15 @@ bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath)
getEdges(eigenEdges);
getVertices(eigenVertices);
vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
label = std::filesystem::path(meshFilePath).stem().string();
label = std::filesystem::path(plyFilename).stem().string();
const bool printDebugInfo = false;
if (printDebugInfo) {
std::cout << meshFilePath << " was loaded successfuly." << std::endl;
std::cout << plyFilename << " was loaded successfuly." << std::endl;
std::cout << "Mesh has " << EN() << " edges." << std::endl;
}
label = std::filesystem::path(meshFilePath).stem().string();
label=std::filesystem::path(plyFilename).stem().string();
return true;
}
@ -248,7 +248,7 @@ bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath)
bool VCGEdgeMesh::loadUsingDefaultLoader(const std::string &plyFilename)
{
Clear();
this->Clear();
// assert(plyFileHasAllRequiredFields(plyFilename));
// Load the ply file
int mask = 0;
@ -337,35 +337,10 @@ bool VCGEdgeMesh::copy(VCGEdgeMesh &mesh) {
return true;
}
void VCGEdgeMesh::set(const std::vector<double> &vertexPositions, const std::vector<int> &edges)
{
Clear();
for (int ei = 0; ei < edges.size(); ei += 2) {
const int vi0 = edges[ei];
const int vi1 = edges[ei + 1];
const int vi0_offset = 3 * vi0;
const int vi1_offset = 3 * vi1;
const CoordType p0(vertexPositions[vi0_offset],
vertexPositions[vi0_offset + 1],
vertexPositions[vi0_offset + 2]);
const CoordType p1(vertexPositions[vi1_offset],
vertexPositions[vi1_offset + 1],
vertexPositions[vi1_offset + 2]);
auto eIt = vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(*this, p0, p1);
CoordType n(0, 0, 1);
eIt->cV(0)->N() = n;
eIt->cV(1)->N() = n;
}
removeDuplicateVertices();
updateEigenEdgeAndVertices();
}
void VCGEdgeMesh::removeDuplicateVertices()
{
vcg::tri::Clean<VCGEdgeMesh>::MergeCloseVertex(*this, 0.000061524);
vcg::tri::Allocator<VCGEdgeMesh>::CompactVertexVector(*this);
vcg::tri::Allocator<VCGEdgeMesh>::CompactEdgeVector(*this);
vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
}
@ -443,12 +418,8 @@ polyscope::CurveNetwork *VCGEdgeMesh::registerForDrawing(
{
PolyscopeInterface::init();
const double drawingRadius = desiredRadius;
updateEigenEdgeAndVertices();
polyscope::CurveNetwork *polyscopeHandle_edgeMesh
= polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges());
// std::cout << "EDGES:" << polyscopeHandle_edgeMesh->nEdges() << std::endl;
assert(polyscopeHandle_edgeMesh->nEdges() == getEigenEdges().rows()
&& polyscopeHandle_edgeMesh->nNodes() == getEigenVertices().rows());
polyscopeHandle_edgeMesh->setEnabled(shouldEnable);
polyscopeHandle_edgeMesh->setRadius(drawingRadius, true);

11
edgemesh.hpp
View File

@ -58,12 +58,9 @@ public:
* */
bool copy(VCGEdgeMesh &mesh);
void set(const std::vector<double> &vertexPositions, const std::vector<int> &edges);
void removeDuplicateVertices();
virtual void deleteDanglingVertices();
virtual void deleteDanglingVertices(
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu);
void deleteDanglingVertices();
void deleteDanglingVertices(vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu);
void getEdges(Eigen::MatrixX3d &edgeStartingPoints, Eigen::MatrixX3d &edgeEndingPoints) const;
@ -73,8 +70,8 @@ public:
// bool loadUsingNanoply(const std::string &plyFilename);
bool load(const std::filesystem::path &meshFilePath) override;
bool save(const std::filesystem::path &meshFilePath = std::filesystem::path()) override;
bool load(const std::string &plyFilename) override;
bool save(const std::string &plyFilename = std::string()) override;
bool createSpanGrid(const size_t squareGridDimension);
bool createSpanGrid(const size_t desiredWidth, const size_t desiredHeight);

2
hexagonremesher.hpp
View File

@ -11,7 +11,7 @@
namespace PolygonalRemeshing {
//std::shared_ptr<VCGPolyMesh>
inline std::shared_ptr<VCGPolyMesh> remeshWithPolygons(VCGTriMesh &surface)
std::shared_ptr<VCGPolyMesh> remeshWithPolygons(VCGTriMesh &surface)
{
// const std::string tileIntoFilePath(
// "/home/iason/Coding/build/PatternTillingReducedModel/RelWithDebInfo/plane.ply");

274
linearsimulationmodel.cpp
View File

@ -1,274 +0,0 @@
#include "linearsimulationmodel.hpp"
//#include "gsl/gsl"
std::vector<fea::Elem> LinearSimulationModel::getFeaElements(
const std::shared_ptr<SimulationMesh> &pMesh)
{
const int numberOfEdges = pMesh->EN();
std::vector<fea::Elem> elements(numberOfEdges);
for (int edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
const SimulationMesh::CoordType &evn0 = pMesh->edge[edgeIndex].cV(0)->cN();
const SimulationMesh::CoordType &evn1 = 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 = pMesh->elements[edgeIndex];
const double E = element.material.youngsModulus;
fea::Props feaProperties(E * element.dimensions.A,
E * element.dimensions.inertia.I3,
E * element.dimensions.inertia.I2,
element.material.G * element.dimensions.inertia.J,
nAverageVector);
const int vi0 = pMesh->getIndex(pMesh->edge[edgeIndex].cV(0));
const int vi1 = pMesh->getIndex(pMesh->edge[edgeIndex].cV(1));
elements[edgeIndex] = fea::Elem(vi0, vi1, feaProperties);
}
return elements;
}
std::vector<fea::Node> LinearSimulationModel::getFeaNodes(
const std::shared_ptr<SimulationMesh> &pMesh)
{
const int numberOfNodes = pMesh->VN();
std::vector<fea::Node> feaNodes(numberOfNodes);
for (int vi = 0; vi < numberOfNodes; vi++) {
const CoordType &p = pMesh->vert[vi].cP();
feaNodes[vi] = fea::Node(p[0], p[1], p[2]);
}
return feaNodes;
}
std::vector<fea::BC> LinearSimulationModel::getFeaFixedNodes(
const std::shared_ptr<SimulationJob> &job)
{
std::vector<fea::BC> boundaryConditions;
// boundaryConditions.reserve(job->constrainedVertices.size() * 6
// + job->nodalForcedDisplacements.size() * 3);
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));
}
}
for (auto forcedDisplacement : job->nodalForcedDisplacements) {
const int vi = forcedDisplacement.first;
for (int dofIndex = 0; dofIndex < 3; dofIndex++) {
boundaryConditions.emplace_back(
fea::BC(vi, static_cast<fea::DOF>(dofIndex), forcedDisplacement.second[dofIndex]));
}
}
return boundaryConditions;
}
std::vector<fea::Force> LinearSimulationModel::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;
}
fea::Force f(nodalForce.first, dofIndex, nodalForce.second[dofIndex]);
nodalForces.emplace_back(f);
}
}
return nodalForces;
}
SimulationResults LinearSimulationModel::getResults(
const fea::Summary &feaSummary, const std::shared_ptr<SimulationJob> &simulationJob)
{
SimulationResults results;
results.converged = feaSummary.converged;
if (!results.converged) {
return results;
}
results.executionTime = feaSummary.total_time_in_ms * 1000;
// displacements
results.displacements.resize(feaSummary.num_nodes);
results.rotationalDisplacementQuaternion.resize(feaSummary.num_nodes);
for (int vi = 0; vi < feaSummary.num_nodes; vi++) {
results.displacements[vi] = Vector6d(feaSummary.nodal_displacements[vi]);
const Vector6d &nodalDisplacement = results.displacements[vi];
Eigen::Quaternion<double> q_nx;
q_nx = Eigen::AngleAxis<double>(nodalDisplacement[3], Eigen::Vector3d(1, 0, 0));
Eigen::Quaternion<double> q_ny;
q_ny = Eigen::AngleAxis<double>(nodalDisplacement[4], Eigen::Vector3d(0, 1, 0));
Eigen::Quaternion<double> q_nz;
q_nz = Eigen::AngleAxis<double>(nodalDisplacement[5], Eigen::Vector3d(0, 0, 1));
results.rotationalDisplacementQuaternion[vi] = q_nx * q_ny * q_nz;
}
results.setLabelPrefix("Linear");
// // 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];
// }
// }
for (int ei = 0; ei < feaSummary.num_elems; ei++) {
const std::vector<double> &elementForces = feaSummary.element_forces[ei];
const Element &element = simulationJob->pMesh->elements[ei];
//Axial
const double elementPotentialEnergy_axial_v0 = std::pow(elementForces[0], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.dimensions.A);
const double elementPotentialEnergy_axial_v1 = std::pow(elementForces[6], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.dimensions.A);
const double elementPotentialEnergy_axial = elementPotentialEnergy_axial_v0
+ elementPotentialEnergy_axial_v1;
//Shear
const double elementPotentialEnergy_shearY_v0 = std::pow(elementForces[1], 2)
* element.initialLength
/ (4 * element.dimensions.A
* element.material.G);
const double elementPotentialEnergy_shearZ_v0 = std::pow(elementForces[2], 2)
* element.initialLength
/ (4 * element.dimensions.A
* element.material.G);
const double elementPotentialEnergy_shearY_v1 = std::pow(elementForces[7], 2)
* element.initialLength
/ (4 * element.dimensions.A
* element.material.G);
const double elementPotentialEnergy_shearZ_v1 = std::pow(elementForces[8], 2)
* element.initialLength
/ (4 * element.dimensions.A
* element.material.G);
const double elementPotentialEnergy_shear = elementPotentialEnergy_shearY_v0
+ elementPotentialEnergy_shearZ_v0
+ elementPotentialEnergy_shearY_v1
+ elementPotentialEnergy_shearZ_v1;
//Bending
const double elementPotentialEnergy_bendingY_v0 = std::pow(elementForces[4], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.dimensions.inertia.I2);
const double elementPotentialEnergy_bendingZ_v0 = std::pow(elementForces[5], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.dimensions.inertia.I3);
const double elementPotentialEnergy_bendingY_v1 = std::pow(elementForces[10], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.dimensions.inertia.I2);
const double elementPotentialEnergy_bendingZ_v1 = std::pow(elementForces[11], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.dimensions.inertia.I3);
const double elementPotentialEnergy_bending = elementPotentialEnergy_bendingY_v0
+ elementPotentialEnergy_bendingZ_v0
+ elementPotentialEnergy_bendingY_v1
+ elementPotentialEnergy_bendingZ_v1;
//Torsion
const double elementPotentialEnergy_torsion_v0 = std::pow(elementForces[3], 2)
* element.initialLength
/ (4 * element.dimensions.inertia.J
* element.material.G);
const double elementPotentialEnergy_torsion_v1 = std::pow(elementForces[9], 2)
* element.initialLength
/ (4 * element.dimensions.inertia.J
* element.material.G);
const double elementPotentialEnergy_torsion = elementPotentialEnergy_torsion_v0
+ elementPotentialEnergy_torsion_v1;
const double elementInternalPotentialEnergy = elementPotentialEnergy_axial
+ elementPotentialEnergy_shear
+ elementPotentialEnergy_bending
+ elementPotentialEnergy_torsion;
results.internalPotentialEnergy += elementInternalPotentialEnergy;
}
results.pJob = simulationJob;
return results;
}
SimulationResults LinearSimulationModel::executeSimulation(
const std::shared_ptr<SimulationJob> &pSimulationJob)
{
assert(pSimulationJob->pMesh->VN() != 0);
const bool wasInitializedWithDifferentStructure = pSimulationJob->pMesh->EN()
!= simulator.structure.elems.size()
|| pSimulationJob->pMesh->VN()
!= simulator.structure.nodes.size();
if (!simulator.wasInitialized() || wasInitializedWithDifferentStructure) {
setStructure(pSimulationJob->pMesh);
}
// 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(pSimulationJob);
auto nodalForces = getFeaNodalForces(pSimulationJob);
fea::Summary feaResults = simulator.solve(fixedVertices, nodalForces, ties, equations, opts);
SimulationResults results = getResults(feaResults, pSimulationJob);
results.pJob = pSimulationJob;
return results;
}
void LinearSimulationModel::setStructure(const std::shared_ptr<SimulationMesh> &pMesh)
{
fea::BeamStructure structure(getFeaNodes(pMesh), getFeaElements(pMesh));
simulator.initialize(structure);
}
void LinearSimulationModel::printInfo(const fea::BeamStructure &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;
}
}

262
linearsimulationmodel.hpp
View File

@ -12,21 +12,265 @@ public:
LinearSimulationModel(){
}
static std::vector<fea::Elem> getFeaElements(const std::shared_ptr<SimulationMesh> &pMesh);
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.inertia.I3,
E * element.inertia.I2,
element.G * element.inertia.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);
}
static std::vector<fea::Node> getFeaNodes(const std::shared_ptr<SimulationMesh> &pMesh);
static std::vector<fea::BC> getFeaFixedNodes(const std::shared_ptr<SimulationJob> &job);
return elements;
}
static std::vector<fea::Force> getFeaNodalForces(const std::shared_ptr<SimulationJob> &job);
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
// + job->nodalForcedDisplacements.size() * 3);
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));
}
}
for (auto forcedDisplacement : job->nodalForcedDisplacements) {
const int vi = forcedDisplacement.first;
for (int dofIndex = 0; dofIndex < 3; dofIndex++) {
boundaryConditions.emplace_back(fea::BC(vi,
static_cast<fea::DOF>(dofIndex),
forcedDisplacement.second[dofIndex]));
}
}
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;
}
fea::Force f(nodalForce.first, dofIndex, nodalForce.second[dofIndex]);
nodalForces.emplace_back(f);
}
}
return nodalForces;
}
static SimulationResults getResults(const fea::Summary &feaSummary,
const std::shared_ptr<SimulationJob> &simulationJob);
const std::shared_ptr<SimulationJob> &simulationJob)
{
SimulationResults results;
results.converged = feaSummary.converged;
if (!results.converged) {
return results;
}
SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &simulationJob);
void setStructure(const std::shared_ptr<SimulationMesh> &pMesh);
results.executionTime = feaSummary.total_time_in_ms * 1000;
// displacements
results.displacements.resize(feaSummary.num_nodes);
results.rotationalDisplacementQuaternion.resize(feaSummary.num_nodes);
for (int vi = 0; vi < feaSummary.num_nodes; vi++) {
results.displacements[vi] = Vector6d(feaSummary.nodal_displacements[vi]);
const Vector6d &nodalDisplacement = results.displacements[vi];
Eigen::Quaternion<double> q_nx;
q_nx = Eigen::AngleAxis<double>(nodalDisplacement[3], Eigen::Vector3d(1, 0, 0));
Eigen::Quaternion<double> q_ny;
q_ny = Eigen::AngleAxis<double>(nodalDisplacement[4], Eigen::Vector3d(0, 1, 0));
Eigen::Quaternion<double> q_nz;
q_nz = Eigen::AngleAxis<double>(nodalDisplacement[5], Eigen::Vector3d(0, 0, 1));
results.rotationalDisplacementQuaternion[vi] = q_nx * q_ny * q_nz;
// results.rotationalDisplacementQuaternion[vi] = q_nz * q_ny * q_nx;
// results.rotationalDisplacementQuaternion[vi] = q_nz * q_nx * q_ny;
}
results.setLabelPrefix("Linear");
// // 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];
// }
// }
for (int ei = 0; ei < feaSummary.num_elems; ei++) {
const std::vector<double> &elementForces = feaSummary.element_forces[ei];
const Element &element = simulationJob->pMesh->elements[ei];
//Axial
const double elementPotentialEnergy_axial_v0 = std::pow(elementForces[0], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.A);
const double elementPotentialEnergy_axial_v1 = std::pow(elementForces[6], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.A);
const double elementPotentialEnergy_axial = elementPotentialEnergy_axial_v0
+ elementPotentialEnergy_axial_v1;
//Shear
const double elementPotentialEnergy_shearY_v0 = std::pow(elementForces[1], 2)
* element.initialLength
/ (4 * element.A * element.G);
const double elementPotentialEnergy_shearZ_v0 = std::pow(elementForces[2], 2)
* element.initialLength
/ (4 * element.A * element.G);
const double elementPotentialEnergy_shearY_v1 = std::pow(elementForces[7], 2)
* element.initialLength
/ (4 * element.A * element.G);
const double elementPotentialEnergy_shearZ_v1 = std::pow(elementForces[8], 2)
* element.initialLength
/ (4 * element.A * element.G);
const double elementPotentialEnergy_shear = elementPotentialEnergy_shearY_v0
+ elementPotentialEnergy_shearZ_v0
+ elementPotentialEnergy_shearY_v1
+ elementPotentialEnergy_shearZ_v1;
//Bending
const double elementPotentialEnergy_bendingY_v0 = std::pow(elementForces[4], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.inertia.I2);
const double elementPotentialEnergy_bendingZ_v0 = std::pow(elementForces[5], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.inertia.I3);
const double elementPotentialEnergy_bendingY_v1 = std::pow(elementForces[10], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.inertia.I2);
const double elementPotentialEnergy_bendingZ_v1 = std::pow(elementForces[11], 2)
* element.initialLength
/ (4 * element.material.youngsModulus
* element.inertia.I3);
const double elementPotentialEnergy_bending = elementPotentialEnergy_bendingY_v0
+ elementPotentialEnergy_bendingZ_v0
+ elementPotentialEnergy_bendingY_v1
+ elementPotentialEnergy_bendingZ_v1;
//Torsion
const double elementPotentialEnergy_torsion_v0 = std::pow(elementForces[3], 2)
* element.initialLength
/ (4 * element.inertia.J * element.G);
const double elementPotentialEnergy_torsion_v1 = std::pow(elementForces[9], 2)
* element.initialLength
/ (4 * element.inertia.J * element.G);
const double elementPotentialEnergy_torsion = elementPotentialEnergy_torsion_v0
+ elementPotentialEnergy_torsion_v1;
const double elementInternalPotentialEnergy = elementPotentialEnergy_axial
+ elementPotentialEnergy_shear
+ elementPotentialEnergy_bending
+ elementPotentialEnergy_torsion;
results.internalPotentialEnergy += elementInternalPotentialEnergy;
}
results.pJob = simulationJob;
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, simulationJob);
results.pJob = simulationJob;
return results;
}
// SimulationResults getResults() const;
// void setResultsNodalDisplacementCSVFilepath(const std::string
// &outputPath); void setResultsElementalForcesCSVFilepath(const std::string
// &outputPath);
private:
fea::ThreedBeamFEA simulator;
static void printInfo(const fea::BeamStructure &job);
// 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

6
mesh.hpp
View File

@ -1,7 +1,7 @@
#ifndef MESH_HPP
#define MESH_HPP
#include <filesystem>
#include <Eigen/Core>
#include <string>
class Mesh
@ -11,8 +11,8 @@ protected:
public:
virtual ~Mesh() = default;
virtual bool load(const std::filesystem::path &meshFilePath) { return false; }
virtual bool save(const std::filesystem::path &meshFilePath) { return false; }
virtual bool load(const std::string &filePath) { return false; }
virtual bool save(const std::string &filePath) { return false; }
std::string getLabel() const;
void setLabel(const std::string &newLabel);

65
polymesh.hpp
View File

@ -5,7 +5,6 @@
#include <filesystem>
#include <wrap/io_trimesh/export_obj.h>
#include <wrap/io_trimesh/export_ply.h>
#include <wrap/io_trimesh/import.h>
#ifdef POLYSCOPE_DEFINED
#include <polyscope/surface_mesh.h>
@ -65,25 +64,14 @@ class VCGPolyMesh
public Mesh
{
public:
virtual bool load(const std::filesystem::__cxx11::path &meshFilePath) override
{
int mask;
vcg::tri::io::Importer<VCGPolyMesh>::LoadMask(meshFilePath.c_str(), mask);
int error = vcg::tri::io::Importer<VCGPolyMesh>::Open(*this, meshFilePath.c_str(), mask);
if (error != 0) {
std::cerr << "Could not load polygonal mesh:" << meshFilePath << std::endl;
return false;
}
vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
vcg::tri::UpdateNormal<VCGPolyMesh>::PerVertexNormalized(*this);
vcg::tri::Clean<VCGPolyMesh>::RemoveUnreferencedVertex(*this);
//finally remove valence 1 vertices on the border
// vcg::PolygonalAlgorithm<PolyMeshType>::RemoveValence2Vertices(dual);
return true;
}
// virtual bool load(const std::string &filename) override {
// int mask;
// vcg::tri::io::Importer<VCGPolyMesh>::LoadMask(filename.c_str(), mask);
// int error = vcg::tri::io::Importer<VCGPolyMesh>::Open(
// *this, filename.c_str(), mask);
// if (error != 0) {
// return false;
// }
// // vcg::tri::io::ImporterOBJ<VCGPolyMesh>::Open();
// // unsigned int mask = 0;
// // mask |= nanoply::NanoPlyWrapper<VCGPolyMesh>::IO_VERTCOORD;
@ -96,6 +84,11 @@ public:
// // 0) {
// // std::cout << "Could not load tri mesh" << std::endl;
// // }
// vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
// vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
// // vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
// vcg::tri::UpdateNormal<VCGPolyMesh>::PerVertexNormalized(*this);
// return true;
// }
Eigen::MatrixX3d getVertices() const
{
@ -120,34 +113,6 @@ public:
return faces;
}
// bool load(const std::filesystem::__cxx11::path &meshFilePath)
// {
// const std::string extension = ".ply";
// std::filesystem::path filePath = meshFilePath;
// assert(std::filesystem::path(filePath).extension().string() == extension);
// unsigned int mask = 0;
// mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
// mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
// mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
// mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
// if (vcg::tri::io::Importer<VCGPolyMesh>::Open(*this, filePath.c_str()) != 0) {
// return false;
// }
// label = meshFilePath.filename();
// return true;
// }
bool save(const std::filesystem::__cxx11::path &meshFilePath = std::filesystem::path())
{
if (meshFilePath.extension() == ".obj") {
return saveOBJ(meshFilePath);
} else if (meshFilePath.extension() == ".ply") {
return savePLY(meshFilePath);
}
return false;
}
bool saveOBJ(const std::filesystem::path &objFilePath = std::filesystem::path())
{
const std::string extension = ".obj";
@ -163,7 +128,7 @@ public:
mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
if (vcg::tri::io::ExporterOBJ<VCGPolyMesh>::Save(*this, filePath.string().c_str(), mask) != 0) {
if (vcg::tri::io::ExporterOBJ<VCGPolyMesh>::Save(*this, filePath.c_str(), mask) != 0) {
return false;
}
return true;
@ -184,7 +149,7 @@ public:
mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
if (vcg::tri::io::ExporterPLY<VCGPolyMesh>::Save(*this, filePath.string().c_str(), mask, false) != 0) {
if (vcg::tri::io::ExporterPLY<VCGPolyMesh>::Save(*this, filePath.c_str(), mask, false) != 0) {
return false;
}
return true;

317
reducedmodelevaluator.cpp
View File

@ -1,317 +0,0 @@
#include "reducedmodelevaluator.hpp"
#include "hexagonremesher.hpp"
#include "trianglepatterngeometry.hpp"
#include <filesystem>
using FullPatternVertexIndex = VertexIndex;
using ReducedPatternVertexIndex = VertexIndex;
ReducedModelEvaluator::ReducedModelEvaluator()
{
}
std::vector<double> ReducedModelEvaluator::evaluateReducedModel(
ReducedModelOptimization::Results &optimizationResults)
{
// std::shared_ptr<VCGPolyMesh> pTileIntoSurface = std::make_shared<VCGPolyMesh>();
// std::filesystem::path tileIntoSurfaceFilePath{
// "/home/iason/Coding/Projects/Approximating shapes with flat "
// "patterns/ReducedModelOptimization/TestSet/TileIntoSurface/plane_34Polygons.ply"};
// assert(std::filesystem::exists(tileIntoSurfaceFilePath));
// pTileIntoSurface->load(tileIntoSurfaceFilePath);
//Set required file paths
const std::filesystem::path tileInto_triMesh_filename
= "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/"
"instantMeshes_plane_34.ply";
const std::filesystem::path reducedPatternFilePath
= "/home/iason/Coding/Projects/Approximating shapes with flat "
"patterns/ReducedModelOptimization/TestSet/ReducedPatterns/single_reduced.ply";
const std::filesystem::path intermediateResultsDirectoryPath
= "/home/iason/Coding/build/ReducedModelOptimization/IntermediateResults";
// const std::filesystem::path drmSettingsFilePath
// = "/home/iason/Coding/Projects/Approximating shapes with flat "
// "patterns/ReducedModelOptimization/DefaultSettings/DRMSettings/"
// "defaultDRMSimulationSettings.json";
DRMSimulationModel::Settings drmSimulationSettings;
// drmSimulationSettings.load(drmSettingsFilePath);
drmSimulationSettings.linearGuessForceScaleFactor = 0.7;
drmSimulationSettings.debugModeStep = 10000;
drmSimulationSettings.beVerbose = true;
constexpr bool shouldRerunFullPatternSimulation = false;
const std::vector<std::string> scenariosTestSetLabels{"22Hex_randomBending0",
"22Hex_randomBending1",
"22Hex_randomBending2",
"22Hex_randomBending3",
"22Hex_randomBending4",
"22Hex_randomBending5",
"22Hex_randomBending6",
"22Hex_randomBending7",
"22Hex_randomBending8",
"22Hex_randomBending9",
"22Hex_randomBending10",
"22Hex_randomBending11",
"22Hex_randomBending12",
"22Hex_randomBending13",
"22Hex_randomBending14",
"22Hex_randomBending15",
"22Hex_randomBending16",
"22Hex_randomBending17",
"22Hex_randomBending18",
"22Hex_randomBending19",
"22Hex_randomBending20",
"22Hex_bending_0.005N",
"22Hex_bending_0.01N",
"22Hex_bending_0.03N",
"22Hex_bending_0.05N",
"22Hex_pullOppositeVerts_0.05N",
"22Hex_pullOppositeVerts_0.1N",
"22Hex_pullOppositeVerts_0.3N",
"22Hex_shear_2N",
"22Hex_shear_5N",
"22Hex_axial_10N",
"22Hex_axial_20N",
"22Hex_cylinder_0.05N",
"22Hex_cylinder_0.1N",
"22Hex_s_0.05N",
"22Hex_s_0.1N"};
//Load surface
std::shared_ptr<VCGPolyMesh> pTileIntoSurface = [&]() {
VCGTriMesh tileInto_triMesh;
const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filename);
assert(surfaceLoadSuccessfull);
return PolygonalRemeshing::remeshWithPolygons(tileInto_triMesh);
}();
const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResults.baseTriangle.cP(0),
(optimizationResults.baseTriangle.cP(1)
+ optimizationResults.baseTriangle.cP(
2))
/ 2);
pTileIntoSurface->moveToCenter();
const double scaleFactor = optimizedBaseTriangleHeight
/ pTileIntoSurface->getAverageFaceRadius();
vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
// tileIntoSurface.registerForDrawing();
// polyscope::show();
//Tile full pattern into surface
std::vector<PatternGeometry> fullPatterns(1);
fullPatterns[0].copy(optimizationResults.baseTriangleFullPattern);
//// Base triangle pattern might contain dangling vertices.Remove those
fullPatterns[0].interfaceNodeIndex = 3;
fullPatterns[0].deleteDanglingVertices();
std::vector<int> perSurfaceFacePatternIndices(pTileIntoSurface->FN(), 0);
std::vector<std::vector<size_t>> perPatternIndexTiledFullPatternEdgeIndices;
std::vector<size_t> tileIntoEdgeToTiledFullVi;
std::shared_ptr<PatternGeometry> pTiledFullPattern
= PatternGeometry::tilePattern(fullPatterns,
{},
*pTileIntoSurface,
perSurfaceFacePatternIndices,
tileIntoEdgeToTiledFullVi,
perPatternIndexTiledFullPatternEdgeIndices);
pTiledFullPattern->setLabel("Tiled_full_patterns");
// pTiledFullPattern->registerForDrawing();
//Tile reduced pattern into surface
PatternGeometry reducedPattern;
reducedPattern.load(reducedPatternFilePath);
reducedPattern.deleteDanglingVertices();
assert(reducedPattern.interfaceNodeIndex == 1);
std::vector<PatternGeometry> reducedPatterns(1);
reducedPatterns[0].copy(reducedPattern);
const auto reducedPatternBaseTriangle = reducedPattern.computeBaseTriangle();
ReducedModelOptimization::Results::applyOptimizationResults_innerHexagon(
optimizationResults, reducedPatternBaseTriangle, reducedPatterns[0]);
std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
std::vector<size_t> tileIntoEdgeToTiledReducedVi;
std::shared_ptr<PatternGeometry> pTiledReducedPattern
= PatternGeometry::tilePattern(reducedPatterns,
{0},
*pTileIntoSurface,
perSurfaceFacePatternIndices,
tileIntoEdgeToTiledReducedVi,
perPatternIndexTiledReducedPatternEdgeIndices);
pTiledReducedPattern->setLabel("Tiled_reduced_patterns");
// pTiledReducedPattern->registerForDrawing();
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
fullToReducedViMap; //of only the common vertices
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
reducedToFullViMap; //of only the common vertices
for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) {
fullToReducedViMap[tileIntoEdgeToTiledFullVi[ei]] = tileIntoEdgeToTiledReducedVi[ei];
}
constructInverseMap(fullToReducedViMap, reducedToFullViMap);
std::vector<size_t> tilledFullPatternInterfaceVi;
tilledFullPatternInterfaceVi.clear();
tilledFullPatternInterfaceVi.resize(fullToReducedViMap.size());
size_t viIndex = 0;
for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
tilledFullPatternInterfaceVi[viIndex++] = fullToReducedPair.first;
}
//Create simulation meshes
////Tessellated full pattern simulation mesh
std::shared_ptr<SimulationMesh> pTiledFullPattern_simulationMesh;
pTiledFullPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledFullPattern);
//NOTE: Those should be derived from the optimization results instead of hardcoding them
pTiledFullPattern_simulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002});
if (optimizationResults.fullPatternYoungsModulus == 0) {
std::cerr << "Full pattern's young modulus not found." << std::endl;
std::terminate();
}
pTiledFullPattern_simulationMesh->setBeamMaterial(0.3,
optimizationResults.fullPatternYoungsModulus);
pTiledFullPattern_simulationMesh->reset();
////Tessellated reduced pattern simulation mesh
std::shared_ptr<SimulationMesh> pTiledReducedPattern_simulationMesh;
pTiledReducedPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledReducedPattern);
const std::vector<size_t> &tiledPatternElementIndicesForReducedPattern
= perPatternIndexTiledReducedPatternEdgeIndices[0];
ReducedModelOptimization::Results::applyOptimizationResults_elements(
optimizationResults, pTiledReducedPattern_simulationMesh);
pTiledReducedPattern_simulationMesh->reset();
std::vector<double> distances_drm2reduced(scenariosTestSetLabels.size(), 0);
std::vector<double> distances_normalizedDrm2reduced(scenariosTestSetLabels.size(), 0);
for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
const std::filesystem::path scenariosDirectoryPath
= "/home/iason/Coding/build/PatternTillingReducedModel/Scenarios/";
const std::filesystem::path tiledReducedPatternJobFilePath
= std::filesystem::path(scenariosDirectoryPath)
.append(jobLabel)
.append("SimulationJobs")
.append("Reduced")
.append(SimulationJob::jsonDefaultFileName);
//set jobs
std::shared_ptr<SimulationJob> pJob_tiledReducedPattern;
pJob_tiledReducedPattern = std::make_shared<SimulationJob>(SimulationJob());
pJob_tiledReducedPattern->load(tiledReducedPatternJobFilePath, false);
pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
std::shared_ptr<SimulationJob> pJob_tiledFullPattern;
pJob_tiledFullPattern = std::make_shared<SimulationJob>(SimulationJob());
pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
pJob_tiledReducedPattern->remap(reducedToFullViMap, *pJob_tiledFullPattern);
// pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern->getLabel());
// pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel());
// polyscope::show();
//Save reduced job
const std::filesystem::path tesellatedResultsFolderPath
= std::filesystem::path(intermediateResultsDirectoryPath).append("TessellatedResults");
const std::filesystem::path surfaceFolderPath = std::filesystem::path(
tesellatedResultsFolderPath)
.append(pTileIntoSurface->getLabel());
const std::string scenarioLabel = pJob_tiledFullPattern->getLabel();
const std::filesystem::path scenarioDirectoryPath = std::filesystem::path(surfaceFolderPath)
.append(scenarioLabel);
const std::filesystem::path reducedJobDirectoryPath
= std::filesystem::path(scenarioDirectoryPath).append("ReducedJob");
std::filesystem::create_directories(reducedJobDirectoryPath);
pJob_tiledReducedPattern->save(reducedJobDirectoryPath);
//Run scenario
////Full
const std::string patternLabel = [&]() {
const std::string patternLabel = optimizationResults.baseTriangleFullPattern.getLabel();
const int numberOfOccurences = std::count_if(patternLabel.begin(),
patternLabel.end(),
[](char c) { return c == '#'; });
if (numberOfOccurences == 0) {
return std::to_string(optimizationResults.baseTriangleFullPattern.EN()) + "#"
+ optimizationResults.baseTriangleFullPattern.getLabel();
} else if (numberOfOccurences == 1) {
return optimizationResults.baseTriangleFullPattern.getLabel();
}
}();
const auto fullResultsFolderPath
= std::filesystem::path(scenarioDirectoryPath).append(patternLabel).append("Results");
if (shouldRerunFullPatternSimulation && std::filesystem::exists(fullResultsFolderPath)) {
std::filesystem::remove_all(fullResultsFolderPath);
}
const std::filesystem::path fullPatternJobFolderPath = std::filesystem::path(
scenarioDirectoryPath)
.append(patternLabel)
.append("SimulationJob");
SimulationResults simulationResults_tiledFullPattern_drm;
if (std::filesystem::exists(fullResultsFolderPath)) {
//Load full pattern results
assert(std::filesystem::exists(fullPatternJobFolderPath));
simulationResults_tiledFullPattern_drm.load(fullResultsFolderPath,
fullPatternJobFolderPath);
simulationResults_tiledFullPattern_drm.converged = true;
} else {
//Full
std::cout << "Executing:" << jobLabel << std::endl;
DRMSimulationModel drmSimulationModel;
simulationResults_tiledFullPattern_drm
= drmSimulationModel.executeSimulation(pJob_tiledFullPattern, drmSimulationSettings);
simulationResults_tiledFullPattern_drm.setLabelPrefix("DRM");
}
std::filesystem::create_directories(fullResultsFolderPath);
simulationResults_tiledFullPattern_drm.save(
std::filesystem::path(scenarioDirectoryPath).append(patternLabel));
if (!simulationResults_tiledFullPattern_drm.converged) {
std::cerr << "Full pattern simulation failed." << std::endl;
}
LinearSimulationModel linearSimulationModel;
SimulationResults simulationResults_tiledReducedPattern
= linearSimulationModel.executeSimulation(pJob_tiledReducedPattern);
// simulationResults_tiledReducedPattern.registerForDrawing();
// simulationResults_tiledFullPattern_drm.registerForDrawing();
// polyscope::show();
//measure distance
const double distance_fullDrmToReduced
= simulationResults_tiledFullPattern_drm
.computeDistance(simulationResults_tiledReducedPattern, fullToReducedViMap);
double distance_fullSumOfAllVerts = 0;
for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
distance_fullSumOfAllVerts += simulationResults_tiledFullPattern_drm
.displacements[fullToReducedPair.first]
.getTranslation()
.norm();
}
const double distance_normalizedFullDrmToReduced = distance_fullDrmToReduced
/ distance_fullSumOfAllVerts;
distances_drm2reduced[jobIndex] = distance_fullDrmToReduced;
distances_normalizedDrm2reduced[jobIndex] = distance_normalizedFullDrmToReduced;
}
//#ifndef POLYSCOPE_DEFINED
// return distances_drm2reduced;
//#else
//report distance
// csvFile csv_results("", flse);
csvFile csv_results({}, false);
// csvFile csv_results(std::filesystem::path(dirPath_thisOptimization)
// .append("results.csv")
// .string(),
// false);
csv_results /*<< "Job Label"*/
<< "drm2Reduced"
<< "norm_drm2Reduced";
csv_results << endrow;
double sumOfNormalizedFull2Reduced = 0;
for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
const double &distance_fullDrmToReduced = distances_drm2reduced[jobIndex];
const double &distance_normalizedFullDrmToReduced = distances_normalizedDrm2reduced[jobIndex];
csv_results /*<< jobLabel*/ << distance_fullDrmToReduced
<< distance_normalizedFullDrmToReduced;
csv_results << endrow;
sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
}
std::cout << "Average normalized error per scenario:"
<< sumOfNormalizedFull2Reduced / scenariosTestSetLabels.size() << std::endl;
return distances_normalizedDrm2reduced;
//#endif
}

14
reducedmodelevaluator.hpp
View File

@ -1,14 +0,0 @@
#ifndef REDUCEDMODELEVALUATOR_HPP
#define REDUCEDMODELEVALUATOR_HPP
#include "reducedmodeloptimizer_structs.hpp"
class ReducedModelEvaluator
{
public:
ReducedModelEvaluator();
static std::vector<double> evaluateReducedModel(
ReducedModelOptimization::Results &optimizationResults);
};
#endif // REDUCEDMODELEVALUATOR_HPP

422
reducedmodeloptimizer_structs.hpp
View File

@ -7,7 +7,7 @@
#include "unordered_map"
#include <string>
namespace ReducedModelOptimization {
namespace ReducedPatternOptimization {
enum BaseSimulationScenario {
Axial,
@ -34,9 +34,9 @@ struct Colors
struct xRange
{
std::string label{};
double min{0};
double max{0};
std::string label;
double min;
double max;
inline bool operator<(const xRange &other) { return label < other.label; }
std::string toString() const
@ -60,106 +60,32 @@ struct xRange
{
return label == xrange.label && min == xrange.min && max == xrange.max;
}
std::tuple<std::string, double, double> toTuple() const
{
return std::make_tuple(label, min, max);
}
void set(const std::tuple<std::string, double, double> &inputTuple)
{
if (std::get<1>(inputTuple) > std::get<2>(inputTuple)) {
std::cerr
<< "Invalid xRange tuple. Second argument(min) cant be smaller than the third(max)"
<< std::endl;
std::terminate();
// return;
}
std::tie(label, min, max) = inputTuple;
}
};
enum OptimizationParameterIndex { E, A, I2, I3, J, R, Theta, NumberOfOptimizationVariables };
inline int getParameterIndex(const std::string &s)
{
if ("E" == s) {
return E;
} else if ("A" == s) {
return A;
} else if ("I2" == s) {
return I2;
} else if ("I3" == s) {
return I3;
} else if ("J" == s) {
return J;
} else if ("R" == s || "HexSize" == s) {
return R;
} else if ("Theta" == s || "HexAngle" == s) {
return Theta;
} else {
std::cerr << "Input is not recognized as a valid optimization variable index:" << s
<< std::endl;
return -1;
}
}
struct Settings
{
inline static std::string defaultFilename{"OptimizationSettings.json"};
std::vector<std::vector<OptimizationParameterIndex>> optimizationStrategy = {
{E, A, I2, I3, J, R, Theta}};
// {A, I2, I3, J, R, Theta}};
std::optional<std::vector<double>>
optimizationVariablesGroupsWeights; //TODO:should be removed in the future if not a splitting strategy is used for the optimization
enum NormalizationStrategy { NonNormalized, Epsilon };
inline static std::vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
NormalizationStrategy normalizationStrategy{Epsilon};
std::array<xRange, NumberOfOptimizationVariables> variablesRanges{xRange{"E", 0.01, 100},
xRange{"A", 0.01, 100},
xRange{"I2", 0.01, 100},
xRange{"I3", 0.01, 100},
xRange{"J", 0.01, 100},
xRange{"R", 0.05, 0.95},
xRange{"Theta", -30, 30}};
std::vector<xRange> xRanges;
inline static vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
int numberOfFunctionCalls{100000};
double solverAccuracy{1e-3};
double translationEpsilon{3e-4};
// double translationEpsilon{0};
double angularDistanceEpsilon{vcg::math::ToRad(0.0)};
double solverAccuracy{1e-2};
NormalizationStrategy normalizationStrategy{Epsilon};
double translationNormalizationParameter{0.003};
double rotationNormalizationParameter{3};
bool splitGeometryMaterialOptimization{false};
struct ObjectiveWeights
{
double translational{1.2};
double rotational{0.8};
bool operator==(const ObjectiveWeights &other) const
{
return this->translational == other.translational
&& this->rotational == other.rotational;
}
};
std::array<ObjectiveWeights, NumberOfBaseSimulationScenarios> perBaseScenarioObjectiveWeights;
std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios>
convertObjectiveWeightsToPairs() const
{
std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios> objectiveWeightsPairs;
for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
baseScenario++) {
objectiveWeightsPairs[baseScenario]
= std::make_pair(perBaseScenarioObjectiveWeights[baseScenario].translational,
perBaseScenarioObjectiveWeights[baseScenario].rotational);
}
return objectiveWeightsPairs;
}
double translational{1};
double rotational{1};
} objectiveWeights;
struct JsonKeys
{
inline static std::string OptimizationStrategy{"OptimizationStrategy"};
inline static std::string OptimizationStrategyGroupWeights{
"OptimizationStrategyGroupWeights"};
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 ObjectiveWeights{"ObjectiveWeight"};
inline static std::string TranslationalObjectiveWeight{"TransObjWeight"};
};
void save(const std::filesystem::path &saveToPath)
@ -167,52 +93,21 @@ struct Settings
assert(std::filesystem::is_directory(saveToPath));
nlohmann::json json;
json[GET_VARIABLE_NAME(optimizationStrategy)] = optimizationStrategy;
if (optimizationVariablesGroupsWeights.has_value()) {
json[GET_VARIABLE_NAME(ptimizationStrategyGroupWeights)]
= optimizationVariablesGroupsWeights.value();
}
//write x ranges
const std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
xRangesAsTuples = [=]() {
std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
xRangesAsTuples;
for (int optimizationParameterIndex = E;
optimizationParameterIndex != NumberOfOptimizationVariables;
optimizationParameterIndex++) {
xRangesAsTuples[optimizationParameterIndex]
= variablesRanges[optimizationParameterIndex].toTuple();
for (size_t xRangeIndex = 0; xRangeIndex < xRanges.size(); xRangeIndex++) {
const auto &xRange = xRanges[xRangeIndex];
json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)]
= xRange.toString();
}
return xRangesAsTuples;
}();
json[JsonKeys::OptimizationVariables] = xRangesAsTuples;
// for (size_t xRangeIndex = 0; xRangeIndex < variablesRanges.size(); xRangeIndex++) {
// const auto &xRange = variablesRanges[xRangeIndex];
// json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)]
// = xRange.toString();
// }
json[GET_VARIABLE_NAME(numberOfFunctionCalls)] = numberOfFunctionCalls;
json[GET_VARIABLE_NAME(solverAccuracy)] = solverAccuracy;
//Objective weights
std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios> objectiveWeightsPairs;
std::transform(perBaseScenarioObjectiveWeights.begin(),
perBaseScenarioObjectiveWeights.end(),
objectiveWeightsPairs.begin(),
[](const ObjectiveWeights &objectiveWeights) {
return std::make_pair(objectiveWeights.translational,
objectiveWeights.rotational);
});
json[JsonKeys::ObjectiveWeights] = objectiveWeightsPairs;
json[GET_VARIABLE_NAME(translationEpsilon)] = translationEpsilon;
json[GET_VARIABLE_NAME(angularDistanceEpsilon)] = vcg::math::ToDeg(
angularDistanceEpsilon);
json[JsonKeys::NumberOfFunctionCalls] = numberOfFunctionCalls;
json[JsonKeys::SolverAccuracy] = solverAccuracy;
json[JsonKeys::TranslationalObjectiveWeight] = objectiveWeights.translational;
std::filesystem::path jsonFilePath(
std::filesystem::path(saveToPath).append(defaultFilename));
std::filesystem::path(saveToPath).append(JsonKeys::filename));
std::ofstream jsonFile(jsonFilePath.string());
jsonFile << json;
jsonFile.close();
}
bool load(const std::filesystem::path &loadFromPath)
@ -221,170 +116,95 @@ struct Settings
//Load optimal X
nlohmann::json json;
std::filesystem::path jsonFilepath(
std::filesystem::path(loadFromPath).append(defaultFilename));
std::filesystem::path(loadFromPath).append(JsonKeys::filename));
if (!std::filesystem::exists(jsonFilepath)) {
std::cerr << "Optimization settings could not be loaded because input filepath does "
"not exist:"
<< jsonFilepath << std::endl;
return false;
}
std::ifstream ifs(jsonFilepath);
ifs >> json;
if (json.contains(GET_VARIABLE_NAME(optimizationStrategy))) {
optimizationStrategy
= std::vector<std::vector<ReducedModelOptimization::OptimizationParameterIndex>>(
(json.at(GET_VARIABLE_NAME(optimizationStrategy))));
}
if (json.contains(GET_VARIABLE_NAME(optimizationStrategyGroupWeights))) {
optimizationVariablesGroupsWeights = std::vector<double>(
json[GET_VARIABLE_NAME(optimizationStrategyGroupWeights)]);
}
//read x ranges
if (json.contains(JsonKeys::OptimizationVariables)) {
const std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
xRangesAsTuples = json.at(JsonKeys::OptimizationVariables);
for (const auto &rangeTuple : xRangesAsTuples) {
variablesRanges[getParameterIndex(std::get<0>(rangeTuple))].set(rangeTuple);
}
} else { //NOTE:legacy compatibility
size_t xRangeIndex = 0;
while (true) {
const std::string jsonXRangeKey = JsonKeys::OptimizationVariables + "_"
+ std::to_string(xRangeIndex++);
+ std::to_string(xRangeIndex);
if (!json.contains(jsonXRangeKey)) {
break;
}
xRange x;
x.fromString(json.at(jsonXRangeKey));
variablesRanges[getParameterIndex(x.label)] = x;
}
xRanges.push_back(x);
xRangeIndex++;
}
if (json.contains(GET_VARIABLE_NAME(numberOfFunctionCalls))) {
numberOfFunctionCalls = json.at(GET_VARIABLE_NAME(numberOfFunctionCalls));
}
if (json.contains(GET_VARIABLE_NAME(solverAccuracy))) {
solverAccuracy = json.at(GET_VARIABLE_NAME(solverAccuracy));
}
//Objective weights
if (json.contains(JsonKeys::ObjectiveWeights)) {
std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios>
objectiveWeightsPairs = json.at(JsonKeys::ObjectiveWeights);
std::transform(objectiveWeightsPairs.begin(),
objectiveWeightsPairs.end(),
perBaseScenarioObjectiveWeights.begin(),
[](const std::pair<double, double> &objectiveWeightsPair) {
return ObjectiveWeights{objectiveWeightsPair.first,
objectiveWeightsPair.second};
});
}
if (json.contains(GET_VARIABLE_NAME(translationalNormalizationEpsilon))) {
translationEpsilon
= json[GET_VARIABLE_NAME(translationalNormalizationEpsilon)];
}
if (json.contains(GET_VARIABLE_NAME(angularDistanceEpsilon))) {
angularDistanceEpsilon = vcg::math::ToRad(
static_cast<double>(json[GET_VARIABLE_NAME(angularDistanceEpsilon)]));
}
// perBaseScenarioObjectiveWeights = json.at(JsonKeys::ObjectiveWeights);
// objectiveWeights.translational = json.at(JsonKeys::ObjectiveWeights);
// objectiveWeights.rotational = 2 - objectiveWeights.translational;
numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls);
solverAccuracy = json.at(JsonKeys::SolverAccuracy);
objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight);
objectiveWeights.rotational = 2 - objectiveWeights.translational;
return true;
}
std::string toString() const
{
std::string settingsString;
if (!variablesRanges.empty()) {
if (!xRanges.empty()) {
std::string xRangesString;
for (const xRange &range : variablesRanges) {
for (const xRange &range : xRanges) {
xRangesString += range.toString() + " ";
}
settingsString += xRangesString;
}
settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls)
+ " Accuracy=" + std::to_string(solverAccuracy);
for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
baseScenario++) {
+" W_t=" + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].translational)
+ " W_r="
+ std::to_string(perBaseScenarioObjectiveWeights[baseScenario].rotational);
}
+ " Accuracy=" + std::to_string(solverAccuracy)
+ " TransWeight=" + std::to_string(objectiveWeights.translational)
+ " RotWeight=" + std::to_string(objectiveWeights.rotational);
return settingsString;
}
void writeHeaderTo(csvFile &os) const
{
if (!variablesRanges.empty()) {
for (const xRange &range : variablesRanges) {
if (!xRanges.empty()) {
for (const xRange &range : xRanges) {
os << range.label + " max";
os << range.label + " min";
}
}
os << "Function Calls";
os << "Solution Accuracy";
os << "Normalization trans epsilon";
os << "Normalization rot epsilon(deg)";
os << JsonKeys::ObjectiveWeights;
os << "Optimization parameters";
os << "Normalization strategy";
os << "Trans weight";
os << "Rot weight";
os << "Splitted geo from mat opt";
// os << std::endl;
}
void writeSettingsTo(csvFile &os) const
{
if (!variablesRanges.empty()) {
for (const xRange &range : variablesRanges) {
if (!xRanges.empty()) {
for (const xRange &range : xRanges) {
os << range.max;
os << range.min;
}
}
os << numberOfFunctionCalls;
os << solverAccuracy;
os << std::to_string(translationEpsilon);
os << std::to_string(vcg::math::ToDeg(angularDistanceEpsilon));
std::string objectiveWeightsString;
objectiveWeightsString += "{";
for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
baseScenario++) {
objectiveWeightsString
+= "{" + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].translational)
+ "," + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].rotational)
+ "}";
}
objectiveWeightsString += "}";
os << objectiveWeightsString;
//export optimization parameters
std::vector<std::vector<int>> vv;
for (const std::vector<OptimizationParameterIndex> &v : optimizationStrategy) {
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);
os << normalizationStrategyStrings[normalizationStrategy] + "_"
+ std::to_string(translationNormalizationParameter);
os << objectiveWeights.translational;
os << objectiveWeights.rotational;
os << (splitGeometryMaterialOptimization == true ? "true" : "false");
}
};
inline bool operator==(const Settings &settings1, const Settings &settings2) noexcept
{
const bool haveTheSameObjectiveWeights
= std::mismatch(settings1.perBaseScenarioObjectiveWeights.begin(),
settings1.perBaseScenarioObjectiveWeights.end(),
settings2.perBaseScenarioObjectiveWeights.begin())
.first
!= settings1.perBaseScenarioObjectiveWeights.end();
return settings1.numberOfFunctionCalls == settings2.numberOfFunctionCalls
&& settings1.variablesRanges == settings2.variablesRanges
&& settings1.solverAccuracy == settings2.solverAccuracy && haveTheSameObjectiveWeights
&& settings1.translationEpsilon
== settings2.translationEpsilon;
&& settings1.xRanges == settings2.xRanges
&& settings1.solverAccuracy == settings2.solverAccuracy
&& settings1.objectiveWeights.translational == settings2.objectiveWeights.translational
&& settings1.translationNormalizationParameter
== settings2.translationNormalizationParameter;
}
inline void updateMeshWithOptimalVariables(const std::vector<double> &x, SimulationMesh &m)
@ -402,9 +222,9 @@ struct Settings
Element &e = m.elements[ei];
e.setDimensions(CrossSectionType(beamWidth, beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
e.dimensions.inertia.J = J;
e.dimensions.inertia.I2 = I2;
e.dimensions.inertia.I3 = I3;
e.inertia.J = J;
e.inertia.I2 = I2;
e.inertia.I3 = I3;
}
CoordType center_barycentric(1, 0, 0);
@ -451,7 +271,6 @@ struct Settings
std::vector<std::pair<std::string, double>> optimalXNameValuePairs;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
double fullPatternYoungsModulus{0};
PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern
vcg::Triangle3<double> baseTriangle;
@ -467,7 +286,6 @@ struct Settings
std::vector<double> perSimulationScenario_translational;
std::vector<double> perSimulationScenario_rotational;
std::vector<double> perSimulationScenario_total;
std::vector<double> perSimulationScenario_total_unweighted;
} objectiveValue;
std::vector<double> perScenario_fullPatternPotentialEnergy;
std::vector<double> objectiveValueHistory;
@ -476,8 +294,8 @@ struct Settings
struct CSVExportingSettings
{
bool exportPE{false};
bool exportIterationOfMinima{false};
bool exportRawObjectiveValue{true};
bool exportIterationOfMinima{true};
bool exportRawObjectiveValue{false};
CSVExportingSettings() {}
};
@ -491,10 +309,9 @@ struct Settings
inline static std::string FullPatternLabel{"FullPatternLabel"};
inline static std::string Settings{"OptimizationSettings"};
inline static std::string FullPatternPotentialEnergies{"PE"};
inline static std::string fullPatternYoungsModulus{"youngsModulus"};
};
void save(const std::string &saveToPath, const bool shouldExportDebugFiles = false)
void save(const string &saveToPath, const bool shouldExportDebugFiles = false)
{
//clear directory
if (std::filesystem::exists(saveToPath)) {
@ -545,7 +362,6 @@ struct Settings
= perScenario_fullPatternPotentialEnergy[simulationScenarioIndex];
}
json_optimizationResults[JsonKeys::FullPatternPotentialEnergies] = fullPatternPE;
json_optimizationResults[JsonKeys::fullPatternYoungsModulus] = fullPatternYoungsModulus;
////Save to json file
std::filesystem::path jsonFilePath(
std::filesystem::path(saveToPath).append(JsonKeys::filename));
@ -567,8 +383,7 @@ struct Settings
= fullPatternSimulationJobs[simulationScenarioIndex];
std::filesystem::path simulationJobFolderPath(
std::filesystem::path(simulationJobsPath)
.append(std::to_string(simulationScenarioIndex) + "_"
+ pFullPatternSimulationJob->label));
.append(pFullPatternSimulationJob->label));
std::filesystem::create_directories(simulationJobFolderPath);
const auto fullPatternDirectoryPath
= std::filesystem::path(simulationJobFolderPath).append("Full");
@ -598,13 +413,12 @@ struct Settings
#endif
}
bool load(const std::string &loadFromPath, const bool &shouldLoadDebugFiles = false)
bool load(const string &loadFromPath, const bool &shouldLoadDebugFiles = false)
{
assert(std::filesystem::is_directory(loadFromPath));
std::filesystem::path jsonFilepath(
std::filesystem::path(loadFromPath).append(JsonKeys::filename));
if (!std::filesystem::exists(jsonFilepath)) {
std::cerr << "Input path does not exist:" << loadFromPath << std::endl;
return false;
}
//Load optimal X
@ -651,83 +465,61 @@ struct Settings
baseTriangle.P2(0) = CoordType(baseTriangleVertices[6],
baseTriangleVertices[7],
baseTriangleVertices[8]);
if (json_optimizationResults.contains(JsonKeys::fullPatternYoungsModulus)) {
fullPatternYoungsModulus = json_optimizationResults.at(
JsonKeys::fullPatternYoungsModulus);
} else {
fullPatternYoungsModulus = 1 * 1e9;
}
if (shouldLoadDebugFiles) {
const std::filesystem::path simulationJobsFolderPath(
std::filesystem::path(loadFromPath).append("SimulationJobs"));
std::vector<std::filesystem::path> sortedByName;
for (auto &entry : std::filesystem::directory_iterator(simulationJobsFolderPath))
sortedByName.push_back(entry.path());
std::sort(sortedByName.begin(), sortedByName.end(), &Utilities::compareNat);
for (const auto &simulationScenarioPath : sortedByName) {
for (const auto &directoryEntry :
filesystem::directory_iterator(simulationJobsFolderPath)) {
const auto simulationScenarioPath = directoryEntry.path();
if (!std::filesystem::is_directory(simulationScenarioPath)) {
continue;
}
// Load full pattern files
for (const auto &fileEntry : std::filesystem::directory_iterator(
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());
job.pMesh->setBeamMaterial(0.3, fullPatternYoungsModulus);
fullPatternSimulationJobs.push_back(std::make_shared<SimulationJob>(job));
}
}
const auto fullJobFilepath = Utilities::getFilepathWithExtension(
std::filesystem::path(simulationScenarioPath).append("Full"), ".json");
SimulationJob fullJob;
fullJob.load(fullJobFilepath.string());
fullJob.pMesh->setBeamMaterial(0.3, fullPatternYoungsModulus);
fullPatternSimulationJobs.push_back(std::make_shared<SimulationJob>(fullJob));
const auto reducedJobFilepath = Utilities::getFilepathWithExtension(
std::filesystem::path(simulationScenarioPath).append("Reduced"), ".json");
SimulationJob reducedJob;
reducedJob.load(reducedJobFilepath.string());
// Load reduced pattern files and apply optimal parameters
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());
// applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh);
applyOptimizationResults_elements(*this, reducedJob.pMesh);
applyOptimizationResults_elements(*this, job.pMesh);
reducedPatternSimulationJobs.push_back(
std::make_shared<SimulationJob>(reducedJob));
std::make_shared<SimulationJob>(job));
}
}
settings.load(std::filesystem::path(loadFromPath).append(Settings::defaultFilename));
}
}
settings.load(loadFromPath);
return true;
}
template<typename MeshType>
static void applyOptimizationResults_innerHexagon(
const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
const vcg::Triangle3<double> &patternBaseTriangle,
MeshType &reducedPattern)
{
std::unordered_map<std::string, double>
optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(),
reducedPattern_optimizationResults.optimalXNameValuePairs.end());
assert(
(optimalXVariables.contains("R") && optimalXVariables.contains("Theta"))
|| (optimalXVariables.contains("HexSize") && optimalXVariables.contains("HexAngle")));
if (optimalXVariables.contains("HexSize")) {
assert(optimalXVariables.contains("HexSize") && optimalXVariables.contains("HexAngle"));
applyOptimizationResults_innerHexagon(optimalXVariables.at("HexSize"),
optimalXVariables.at("HexAngle"),
patternBaseTriangle,
reducedPattern);
return;
}
applyOptimizationResults_innerHexagon(optimalXVariables.at("R"),
optimalXVariables.at("Theta"),
patternBaseTriangle,
reducedPattern);
}
template<typename MeshType>
@ -768,8 +560,8 @@ struct Settings
}
static void applyOptimizationResults_elements(
const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
const std::shared_ptr<SimulationMesh> &pReducedPattern_simulationMesh)
const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
const shared_ptr<SimulationMesh> &pTiledReducedPattern_simulationMesh)
{
assert(CrossSectionType::name == RectangularBeamDimensions::name);
// Set reduced parameters fron the optimization results
@ -783,8 +575,8 @@ struct Settings
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
CrossSectionType elementDimensions(beamWidth, beamHeight);
for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pReducedPattern_simulationMesh->elements[ei];
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.setDimensions(elementDimensions);
}
}
@ -794,8 +586,8 @@ struct Settings
if (optimalXVariables.contains(ymLabel)) {
const double E = optimalXVariables.at(ymLabel);
const ElementMaterial elementMaterial(poissonsRatio, E);
for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pReducedPattern_simulationMesh->elements[ei];
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.setMaterial(elementMaterial);
}
}
@ -803,30 +595,30 @@ struct Settings
const std::string JLabel = "J";
if (optimalXVariables.contains(JLabel)) {
const double J = optimalXVariables.at(JLabel);
for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pReducedPattern_simulationMesh->elements[ei];
e.dimensions.inertia.J = J;
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.inertia.J = J;
}
}
const std::string I2Label = "I2";
if (optimalXVariables.contains(I2Label)) {
const double I2 = optimalXVariables.at(I2Label);
for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pReducedPattern_simulationMesh->elements[ei];
e.dimensions.inertia.I2 = I2;
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.inertia.I2 = I2;
}
}
const std::string I3Label = "I3";
if (optimalXVariables.contains(I3Label)) {
const double I3 = optimalXVariables.at(I3Label);
for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pReducedPattern_simulationMesh->elements[ei];
e.dimensions.inertia.I3 = I3;
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.inertia.I3 = I3;
}
}
pReducedPattern_simulationMesh->reset();
pTiledReducedPattern_simulationMesh->reset();
}
#if POLYSCOPE_DEFINED
@ -923,7 +715,7 @@ struct Settings
}
}
void writeHeaderTo(csvFile &os) const
void writeHeaderTo(csvFile &os)
{
if (exportSettings.exportRawObjectiveValue) {
os << "Total raw Obj value";
@ -959,20 +751,13 @@ struct Settings
// os << "notes";
}
void writeHeaderTo(std::vector<csvFile *> &vectorOfPointersToOutputStreams) const
{
for (csvFile *outputStream : vectorOfPointersToOutputStreams) {
writeHeaderTo(*outputStream);
}
}
void writeResultsTo(csvFile &os) const
{
if (exportSettings.exportRawObjectiveValue) {
os << objectiveValue.totalRaw;
}
os << objectiveValue.total;
if (exportSettings.exportIterationOfMinima && !objectiveValueHistory_iteration.empty()) {
if (exportSettings.exportIterationOfMinima) {
os << objectiveValueHistory_iteration.back();
}
for (int simulationScenarioIndex = 0;
@ -1002,13 +787,6 @@ struct Settings
// os << notes;
}
void writeResultsTo(std::vector<csvFile *> &vectorOfPointersToOutputStreams) const
{
for (csvFile *&outputStream : vectorOfPointersToOutputStreams) {
writeResultsTo(*outputStream);
}
}
void writeMinimaInfoTo(csvFile &outputCsv)
{
outputCsv << "Iteration";

243
simulation_structs.hpp
View File

@ -1,6 +1,5 @@
#ifndef SIMULATIONSTRUCTS_HPP
#define SIMULATIONSTRUCTS_HPP
#include <fstream>
namespace Eigen {
template <class Matrix>
@ -24,13 +23,13 @@ void read_binary(const std::string &filename, Matrix &matrix) {
in.read((char *)matrix.data(), rows * cols * sizeof(typename Matrix::Scalar));
in.close();
}
//const static IOFormat CSVFormat(StreamPrecision, DontAlignCols, ", ", "\n");
//template<class Matrix>
//void writeToCSV(const std::string &filename, Matrix &matrix)
//{
// ofstream file(filename.c_str());
// file << matrix.format(CSVFormat);
//}
const static IOFormat CSVFormat(StreamPrecision, DontAlignCols, ", ", "\n");
template<class Matrix>
void writeToCSV(const std::string &filename, Matrix &matrix)
{
ofstream file(filename.c_str());
file << matrix.format(CSVFormat);
}
} // namespace Eigen
#include "simulationmesh.hpp"
@ -113,20 +112,18 @@ class SimulationJob {
} jsonLabels;
public:
inline static std::string jsonDefaultFileName = "SimulationJob.json";
std::shared_ptr<SimulationMesh> pMesh;
std::string label{"empty_job"};
std::unordered_map<VertexIndex, std::unordered_set<int>> constrainedVertices;
std::unordered_map<VertexIndex, Vector6d> nodalExternalForces;
std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
SimulationJob(const std::shared_ptr<SimulationMesh> &m,
const std::string &label,
SimulationJob(
const std::shared_ptr<SimulationMesh> &m, const std::string &label,
const std::unordered_map<VertexIndex, std::unordered_set<int>> &cv,
const std::unordered_map<VertexIndex, Vector6d> &ef = {},
const std::unordered_map<VertexIndex, Eigen::Vector3d> &fd = {})
: pMesh(m), label(label), constrainedVertices(cv), nodalExternalForces(ef),
nodalForcedDisplacements(fd)
{}
: pMesh(m), label(label), constrainedVertices(cv),
nodalExternalForces(ef), nodalForcedDisplacements(fd) {}
SimulationJob() {}
SimulationJob(const std::string &jsonFilename) { load(jsonFilename); }
@ -135,32 +132,32 @@ public:
return constrainedVertices.empty() && nodalExternalForces.empty()
&& nodalForcedDisplacements.empty() && pMesh == nullptr;
}
void remap(const std::unordered_map<size_t, size_t> &sourceToDestinationViMap,
SimulationJob &destination_simulationJob) const
void remap(const std::unordered_map<size_t, size_t> &thisToOtherViMap,
SimulationJob &simulationJobMapped) const
{
std::unordered_map<VertexIndex, std::unordered_set<int>> destination_fixedVertices;
for (const auto &source_fixedVertex : this->constrainedVertices) {
destination_fixedVertices[sourceToDestinationViMap.at(source_fixedVertex.first)]
= source_fixedVertex.second;
assert(simulationJobMapped.pMesh->VN() != 0);
std::unordered_map<VertexIndex, std::unordered_set<int>> reducedModelFixedVertices;
for (auto fullModelFixedVertex : this->constrainedVertices) {
reducedModelFixedVertices[thisToOtherViMap.at(fullModelFixedVertex.first)]
= fullModelFixedVertex.second;
}
std::unordered_map<VertexIndex, Vector6d> destination_nodalForces;
for (const auto &source_nodalForces : this->nodalExternalForces) {
destination_nodalForces[sourceToDestinationViMap.at(source_nodalForces.first)]
= source_nodalForces.second;
std::unordered_map<VertexIndex, Vector6d> reducedModelNodalForces;
for (auto fullModelNodalForce : this->nodalExternalForces) {
reducedModelNodalForces[thisToOtherViMap.at(fullModelNodalForce.first)]
= fullModelNodalForce.second;
}
std::unordered_map<VertexIndex, Eigen::Vector3d> destination_forcedDisplacements;
for (const auto &source_forcedDisplacements : this->nodalForcedDisplacements) {
destination_forcedDisplacements[sourceToDestinationViMap.at(
source_forcedDisplacements.first)]
= source_forcedDisplacements.second;
std::unordered_map<VertexIndex, Eigen::Vector3d> reducedNodalForcedDisplacements;
for (auto fullForcedDisplacement : this->nodalForcedDisplacements) {
reducedNodalForcedDisplacements[thisToOtherViMap.at(fullForcedDisplacement.first)]
= fullForcedDisplacement.second;
}
destination_simulationJob.constrainedVertices = destination_fixedVertices;
destination_simulationJob.nodalExternalForces = destination_nodalForces;
destination_simulationJob.label = this->getLabel();
destination_simulationJob.nodalForcedDisplacements = destination_forcedDisplacements;
simulationJobMapped.constrainedVertices = reducedModelFixedVertices;
simulationJobMapped.nodalExternalForces = reducedModelNodalForces;
simulationJobMapped.label = this->getLabel();
simulationJobMapped.nodalForcedDisplacements = reducedNodalForcedDisplacements;
}
SimulationJob getCopy() const {
SimulationJob jobCopy;
@ -379,7 +376,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
return;
}
std::vector<std::array<double, 3>> nodeColors(pMesh->VN());
for (const auto &fixedVertex : constrainedVertices) {
for (auto fixedVertex : constrainedVertices) {
const bool hasRotationalDoFConstrained = fixedVertex.second.contains(3)
|| fixedVertex.second.contains(4)
|| fixedVertex.second.contains(5);
@ -395,8 +392,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
}
}
if (!nodalForcedDisplacements.empty()) {
for (const std::pair<VertexIndex, Eigen::Vector3d> &viDisplPair :
nodalForcedDisplacements) {
for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair : nodalForcedDisplacements) {
const VertexIndex vi = viDisplPair.first;
nodeColors[vi][0] += 1;
nodeColors[vi][0] /= 2;
@ -491,7 +487,6 @@ struct SimulationResults
/*TODO: remove rotationalDisplacementQuaternion since the last three components of the displacments
* vector contains the same info using euler angles
*/
inline const static std::string defaultJsonFilename{"SimulationResults.json"};
bool converged{false};
std::shared_ptr<SimulationJob> pJob;
SimulationHistory history;
@ -543,7 +538,7 @@ struct SimulationResults
const std::filesystem::path outputFolderPath = outputFolder.empty()
? std::filesystem::current_path()
: 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(outputFolderPath).append("SimulationJob");
std::filesystem::create_directories(simulationJobOutputFolderPath);
@ -556,16 +551,6 @@ struct SimulationResults
std::filesystem::create_directories(resultsFolderPath);
Eigen::write_binary(std::filesystem::path(resultsFolderPath).append(filename).string(), m);
nlohmann::json json;
json[GET_VARIABLE_NAME(internalPotentialEnergy)] = internalPotentialEnergy;
std::filesystem::path jsonFilePath(
std::filesystem::path(resultsFolderPath).append(defaultJsonFilename));
std::ofstream jsonFile(jsonFilePath.string());
jsonFile << json;
jsonFile.close();
saveDeformedModel(resultsFolderPath.string());
}
@ -583,41 +568,34 @@ struct SimulationResults
void load(const std::filesystem::path &loadFromPath, const std::filesystem::path &loadJobFrom)
{
//load job
pJob->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json").string());
load(loadFromPath);
//Use the first .eigenBin file for loading the displacements
for (auto const &entry : std::filesystem::recursive_directory_iterator(loadFromPath)) {
if (filesystem::is_regular_file(entry) && entry.path().extension() == ".eigenBin") {
Eigen::MatrixXd displacements_eigen;
Eigen::read_binary(entry.path().string(), displacements_eigen);
displacements = Utilities::fromEigenMatrix(displacements_eigen);
break;
}
void load(const std::filesystem::path &loadFromPath, const std::shared_ptr<SimulationJob> &pJob)
{
this->pJob = pJob;
load(loadFromPath);
}
static double computeDistance(
const SimulationResults &resultsA,
const SimulationResults &resultsB,
const std::unordered_map<VertexIndex, VertexIndex> &resultsAToResultsBViMap)
{
double distance = 0;
for (std::pair<int, int> resultsAToResultsBViPair : resultsAToResultsBViMap) {
const double vertexToVertexDistance
= (resultsA.displacements[resultsAToResultsBViPair.first].getTranslation()
- resultsB.displacements[resultsAToResultsBViPair.second].getTranslation())
.norm();
distance += vertexToVertexDistance;
}
return distance;
}
double computeDistance(const SimulationResults &other,
const std::unordered_map<VertexIndex, VertexIndex> &thisToOtherViMap) const
{
return computeDistance(*this, other, thisToOtherViMap);
rotationalDisplacementQuaternion.resize(pJob->pMesh->VN());
for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
rotationalDisplacementQuaternion[vi] = Eigen::AngleAxisd(displacements[vi][3],
Eigen::Vector3d::UnitX())
* Eigen::AngleAxisd(displacements[vi][4],
Eigen::Vector3d::UnitY())
* Eigen::AngleAxisd(displacements[vi][5],
Eigen::Vector3d::UnitZ());
}
}
#ifdef POLYSCOPE_DEFINED
void unregister() const
{
void unregister() const {
if (!polyscope::hasCurveNetwork(getLabel())) {
std::cerr << "No curve network registered with a name: " << getLabel() << std::endl;
std::cerr << "No curve network registered with a name: " << getLabel()
<< std::endl;
std::cerr << "Nothing to remove." << std::endl;
return;
}
@ -635,13 +613,9 @@ struct SimulationResults
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());
mesh->getEigenVertices(),
mesh->getEigenEdges());
}
polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable);
polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, true);
@ -661,18 +635,13 @@ struct SimulationResults
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{};
std::unordered_set<int> interfaceNodes{3, 7, 11, 15, 19, 23};
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)) {
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);
@ -688,14 +657,14 @@ struct SimulationResults
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);
// }
} 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);
@ -723,10 +692,10 @@ struct SimulationResults
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_initialMesh = polyscope::getCurveNetwork(mesh->getLabel());
// if (!polyscopeHandle_initialMesh) {
// polyscopeHandle_initialMesh = mesh->registerForDrawing();
// }
// auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh
// ->addNodeVectorQuantity("FrameX", framesX_initial);
@ -748,66 +717,24 @@ struct SimulationResults
// /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
// //}
pJob->registerForDrawing(getLabel());
// static bool wasExecuted =false;
// if (!wasExecuted) {
// std::function<void()> callback = [&]() {
// static bool showFrames = shouldShowFrames;
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:
bool 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
for (auto const &entry : std::filesystem::recursive_directory_iterator(loadFromPath)) {
if (std::filesystem::is_regular_file(entry) && entry.path().extension() == ".eigenBin") {
Eigen::MatrixXd displacements_eigen;
Eigen::read_binary(entry.path().string(), displacements_eigen);
displacements = Utilities::fromEigenMatrix(displacements_eigen);
break;
}
}
rotationalDisplacementQuaternion.resize(pJob->pMesh->VN());
for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
rotationalDisplacementQuaternion[vi] = Eigen::AngleAxisd(displacements[vi][3],
Eigen::Vector3d::UnitX())
* Eigen::AngleAxisd(displacements[vi][4],
Eigen::Vector3d::UnitY())
* Eigen::AngleAxisd(displacements[vi][5],
Eigen::Vector3d::UnitZ());
}
const std::filesystem::path jsonFilepath = std::filesystem::path(loadFromPath)
.append(defaultJsonFilename);
if (!std::filesystem::exists(jsonFilepath)) {
std::cerr << "Simulation results could not be loaded because filepath does "
"not exist:"
<< jsonFilepath << std::endl;
return false;
}
std::ifstream ifs(jsonFilepath);
nlohmann::json json;
ifs >> json;
if (json.contains(GET_VARIABLE_NAME(internalPotentialEnergy))) {
internalPotentialEnergy = json.at(GET_VARIABLE_NAME(internalPotentialEnergy));
}
return true;
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

32
simulationhistoryplotter.hpp
View File

@ -1,8 +1,8 @@
#ifndef SIMULATIONHISTORYPLOTTER_HPP
#define SIMULATIONHISTORYPLOTTER_HPP
#include "simulation_structs.hpp"
#include "simulationmesh.hpp"
#include "simulation_structs.hpp"
#include "utilities.hpp"
#include <algorithm>
#include <matplot/matplot.h>
@ -15,8 +15,10 @@ struct SimulationResultsReporter {
void writeStatistics(const SimulationResults &results,
const std::string &reportFolderPath) {
std::ofstream file;
file.open(std::filesystem::path(reportFolderPath).append("results.txt").string());
ofstream file;
file.open(
std::filesystem::path(reportFolderPath).append("results.txt").string());
const size_t numberOfSteps = results.history.numberOfSteps;
file << "Number of steps " << numberOfSteps << "\n";
@ -112,23 +114,35 @@ struct SimulationResultsReporter {
if (YvaluesToPlot.size() < 2) {
return;
}
std::vector<double> colors(YvaluesToPlot.size(), 0.5);
std::vector<double> markerSizes(YvaluesToPlot.size(), 5);
auto x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size());
std::vector<double> colors(x.size(), 0.5);
std::vector<double> markerSizes(x.size(), 5);
if (!markPoints.empty()) {
for (const auto pointIndex : markPoints) {
colors[pointIndex] = 0.9;
markerSizes[pointIndex] = 14;
}
}
std::vector<double> x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size());
createPlot(xLabel, yLabel, x, YvaluesToPlot, markerSizes, colors, saveTo);
// matplot::figure(true);
// matplot::hold(matplot::on);
// ->marker_indices(history.redMarks)
// ->marker_indices(truncatedRedMarks)
// .marker_color("r")
// ->marker_size(1)
// ;
// auto greenMarksY = matplot::transform(history.greenMarks,
// [&](auto x) { return YvaluesToPlot[x]; });
// matplot::scatter(history.greenMarks, greenMarksY)->color("green").marker_size(10);
// matplot::hold(matplot::off);
}
void createPlots(const SimulationHistory &history,
const std::string &reportFolderPath,
const std::string &graphSuffix)
{
const auto graphsFolder = std::filesystem::path(reportFolderPath).append("Graphs");
const std::string &graphSuffix) {
const auto graphsFolder =
std::filesystem::path(reportFolderPath).append("Graphs");
std::filesystem::remove_all(graphsFolder);
std::filesystem::create_directory(graphsFolder.string());

94
simulationmesh.cpp
View File

@ -15,8 +15,8 @@ SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) {
elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(
*this, std::string("Elements"));
nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(*this,
std::string("Nodes"));
nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(
*this, std::string("Nodes"));
initializeNodes();
initializeElements();
}
@ -198,6 +198,7 @@ void SimulationMesh::reset() {
}
void SimulationMesh::initializeElements() {
computeElementalProperties();
for (const EdgeType &e : edge) {
Element &element = elements[e];
element.ei = getIndex(e);
@ -242,6 +243,8 @@ void SimulationMesh::updateElementalLengths() {
const vcg::Segment3<double> s(p0, p1);
const double elementLength = s.Length();
elements[e].length = elementLength;
int i = 0;
i++;
}
}
@ -257,6 +260,7 @@ void SimulationMesh::setBeamCrossSection(
const CrossSectionType &beamDimensions) {
for (size_t ei = 0; ei < EN(); ei++) {
elements[ei].dimensions = beamDimensions;
elements[ei].computeDimensionsProperties(beamDimensions);
elements[ei].updateRigidity();
}
}
@ -284,8 +288,7 @@ std::vector<ElementMaterial> SimulationMesh::getBeamMaterial() {
return beamMaterial;
}
bool SimulationMesh::load(const std::string &plyFilename)
{
bool SimulationMesh::load(const string &plyFilename) {
this->Clear();
// assert(plyFileHasAllRequiredFields(plyFilename));
// Load the ply file
@ -325,7 +328,7 @@ bool SimulationMesh::load(const std::string &plyFilename)
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
// mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_MESHATTRIB;
if (!VCGEdgeMesh::load(plyFilename)) {
if (!load(plyFilename)) {
return false;
}
@ -336,7 +339,6 @@ bool SimulationMesh::load(const std::string &plyFilename)
vcg::tri::UpdateTopology<SimulationMesh>::VertexEdge(*this);
initializeNodes();
initializeElements();
setBeamMaterial(0.3, 1 * 1e9);
updateEigenEdgeAndVertices();
// if (!handleBeamProperties._handle->data.empty()) {
@ -356,7 +358,8 @@ bool SimulationMesh::load(const std::string &plyFilename)
if (normalsAreAbsent) {
CoordType normalVector(0, 0, 1);
std::cout << "Warning: Normals are missing from " << plyFilename
<< ". Added normal vector:" << toString(normalVector) << std::endl;
<< ". Added normal vector:" << toString(normalVector)
<< std::endl;
for (auto &v : vert) {
v.N() = normalVector;
}
@ -395,7 +398,7 @@ bool SimulationMesh::save(const std::string &plyFilename)
// != 1) {
// return false;
// }
if (!VCGEdgeMesh::save(plyFilename)) {
if (!save(plyFilename)) {
return false;
}
return true;
@ -453,59 +456,60 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1,
return angle;
}
//void Element::computeMaterialProperties(const ElementMaterial &material) {
// G = material.youngsModulus / (2 * (1 + material.poissonsRatio));
//}
void Element::computeMaterialProperties(const ElementMaterial &material) {
G = material.youngsModulus / (2 * (1 + material.poissonsRatio));
}
//void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A)
//{
// A = dimensions.b * dimensions.h;
//}
void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A)
{
A = dimensions.b * dimensions.h;
}
//void Element::computeDimensionsProperties(
// const RectangularBeamDimensions &dimensions) {
// assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions));
// computeCrossSectionArea(dimensions, A);
// computeMomentsOfInertia(dimensions, dimensions.inertia);
//}
void Element::computeDimensionsProperties(
const RectangularBeamDimensions &dimensions) {
assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions));
computeCrossSectionArea(dimensions, A);
computeMomentsOfInertia(dimensions, inertia);
}
//void Element::computeDimensionsProperties(
// const CylindricalBeamDimensions &dimensions) {
// assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
// A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2));
// dimensions.inertia.I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64;
// dimensions.inertia.I3 = dimensions.inertia.I2;
// dimensions.inertia.J = dimensions.inertia.I2 + dimensions.inertia.I3;
//}
void Element::computeDimensionsProperties(
const CylindricalBeamDimensions &dimensions) {
assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2));
inertia.I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64;
inertia.I3 = inertia.I2;
inertia.J = inertia.I2 + inertia.I3;
}
void Element::setDimensions(const CrossSectionType &dimensions) {
this->dimensions = dimensions;
assert(this->dimensions.A == dimensions.A);
// computeDimensionsProperties(dimensions);
computeDimensionsProperties(dimensions);
updateRigidity();
}
void Element::setMaterial(const ElementMaterial &material)
{
void Element::setMaterial(const ElementMaterial &material) {
this->material = material;
// computeMaterialProperties(material);
computeMaterialProperties(material);
updateRigidity();
}
double Element::getMass(const double &materialDensity)
{
const double beamVolume = dimensions.A * length;
const double beamVolume = A * length;
return beamVolume * materialDensity;
}
void Element::updateRigidity() {
// assert(initialLength != 0);
rigidity.axial = material.youngsModulus * dimensions.A / initialLength;
// assert(rigidity.axial != 0);
rigidity.torsional = material.G * dimensions.inertia.J / initialLength;
// assert(rigidity.torsional != 0);
rigidity.firstBending = 2 * material.youngsModulus * dimensions.inertia.I2 / initialLength;
// assert(rigidity.firstBending != 0);
rigidity.secondBending = 2 * material.youngsModulus * dimensions.inertia.I3 / initialLength;
// assert(rigidity.secondBending != 0);
void Element::computeMomentsOfInertia(const RectangularBeamDimensions &dimensions,
Element::MomentsOfInertia &inertia)
{
inertia.I2 = dimensions.b * std::pow(dimensions.h, 3) / 12;
inertia.I3 = dimensions.h * std::pow(dimensions.b, 3) / 12;
inertia.J = inertia.I2 + inertia.I3;
}
void Element::updateRigidity() {
rigidity.axial = material.youngsModulus * A / initialLength;
rigidity.torsional = G * inertia.J / initialLength;
rigidity.firstBending = 2 * material.youngsModulus * inertia.I2 / initialLength;
rigidity.secondBending = 2 * material.youngsModulus * inertia.I3 / initialLength;
}

20
simulationmesh.hpp
View File

@ -1,9 +1,9 @@
#ifndef SIMULATIONMESH_HPP
#define SIMULATIONMESH_HPP
#include "Eigen/Dense"
#include "edgemesh.hpp"
#include "trianglepatterngeometry.hpp"
#include <Eigen/Dense>
//extern bool drawGlobal;
struct Element;
@ -34,7 +34,7 @@ public:
std::vector<VCGEdgeMesh::EdgePointer>
getIncidentElements(const VCGEdgeMesh::VertexType &v);
virtual bool load(const std::string &plyFilename);
virtual bool load(const string &plyFilename);
std::vector<CrossSectionType> getBeamDimensions();
std::vector<ElementMaterial> getBeamMaterial();
@ -62,14 +62,26 @@ struct Element {
CrossSectionType dimensions;
ElementMaterial material;
double G{0};
double A{0}; // cross sectional area
void computeMaterialProperties(const ElementMaterial &material);
// void computeDimensionsProperties(const RectangularBeamDimensions &dimensions);
// void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions);
void computeDimensionsProperties(const RectangularBeamDimensions &dimensions);
void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions);
void setDimensions(const CrossSectionType &dimensions);
void setMaterial(const ElementMaterial &material);
double getMass(const double &matrialDensity);
struct MomentsOfInertia
{
double I2{0}; // second moment of inertia
double I3{0}; // third moment of inertia
double J{0}; // torsional constant (polar moment of inertia)
} inertia;
static void computeMomentsOfInertia(const RectangularBeamDimensions &dimensions,
MomentsOfInertia &inertia);
struct LocalFrame {
VectorType t1;
VectorType t2;

35
trianglepatterngeometry.cpp
View File

@ -86,8 +86,8 @@ std::vector<vcg::Point3d> PatternGeometry::getVertices() const {
return verts;
}
PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern)
{
PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern) {
const size_t fanSize = PatternGeometry().getFanSize();
PatternGeometry fan(createFan(pattern));
PatternGeometry tile(fan);
@ -99,15 +99,11 @@ PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern)
vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, R);
}
vcg::Matrix44d T;
// const double centerAngle = 2 * M_PI / fanSize;
// const double triangleHeight = std::sin((M_PI - centerAngle) / 2)
// * PatternGeometry().triangleEdgeSize;
vcg::tri::UpdateBounding<PatternGeometry>::Box(fan);
const double triangleHeight = fan.bbox.DimY() / 2;
const double centerAngle = 2 * M_PI / fanSize;
const double triangleHeight =
std::sin((M_PI - centerAngle) / 2) * PatternGeometry().triangleEdgeSize;
T.SetTranslate(0, -2 * triangleHeight, 0);
vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, T);
// fan.registerForDrawing();
// polyscope::show();
PatternGeometry fanOfFan = createFan(fan);
vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(tile, fanOfFan);
@ -149,9 +145,9 @@ PatternGeometry::PatternGeometry(PatternGeometry &other) {
vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
}
bool PatternGeometry::load(const std::filesystem::__cxx11::path &meshFilePath)
bool PatternGeometry::load(const string &plyFilename)
{
if (!VCGEdgeMesh::load(meshFilePath)) {
if (!VCGEdgeMesh::load(plyFilename)) {
return false;
}
vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
@ -781,22 +777,6 @@ double PatternGeometry::computeBaseTriangleHeight() const
return vcg::Distance(vert[0].cP(), vert[interfaceNodeIndex].cP());
}
void PatternGeometry::deleteDanglingVertices()
{
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> pu;
VCGEdgeMesh::deleteDanglingVertices(pu);
if (!pu.remap.empty()) {
interfaceNodeIndex
= pu.remap[interfaceNodeIndex]; //TODO:Could this be automatically be determined?
}
}
void PatternGeometry::deleteDanglingVertices(
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu)
{
VCGEdgeMesh::deleteDanglingVertices(pu);
}
vcg::Triangle3<double> PatternGeometry::getBaseTriangle() const
{
return baseTriangle;
@ -834,7 +814,6 @@ PatternGeometry::PatternGeometry(
updateEigenEdgeAndVertices();
}
//TODO: refactor such that it takes as an input a single pattern and tiles into the desired faces without requiring the each face will contain a single pattern
std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
std::vector<ConstPatternGeometry> &patterns,
const std::vector<int> &connectToNeighborsVi,

5
trianglepatterngeometry.hpp
View File

@ -46,7 +46,7 @@ private:
* vcglib interface.
* */
PatternGeometry(PatternGeometry &other);
bool load(const std::filesystem::path &meshFilePath) override;
bool load(const std::string &plyFilename) override;
void add(const std::vector<vcg::Point3d> &vertices);
void add(const std::vector<vcg::Point2i> &edges);
void add(const std::vector<vcg::Point3d> &vertices, const std::vector<vcg::Point2i> &edges);
@ -124,9 +124,6 @@ private:
std::vector<std::vector<size_t>> &perPatternIndexTiledPatternEdgeIndex);
std::unordered_set<VertexIndex> getInterfaceNodes(const std::vector<size_t> &numberOfNodesPerSlot);
bool isInterfaceConnected(const std::unordered_set<VertexIndex> &interfaceNodes);
void deleteDanglingVertices() override;
void deleteDanglingVertices(
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu) override;
};
#endif // FLATPATTERNGEOMETRY_HPP

159
utilities.hpp
View File

@ -2,14 +2,12 @@
#define UTILITIES_H
#include <Eigen/Dense>
#include <algorithm>
#include <array>
#include <chrono>
#include <filesystem>
#include <fstream>
#include <iterator>
#include <numeric>
#include <regex>
#include <iterator>
#include <algorithm>
#include <array>
#define GET_VARIABLE_NAME(Variable) (#Variable)
@ -142,110 +140,7 @@ struct Vector6d : public std::array<double, 6> {
};
namespace Utilities {
inline bool compareNat(const std::string &a, const std::string &b)
{
if (a.empty())
return true;
if (b.empty())
return false;
if (std::isdigit(a[0]) && !std::isdigit(b[0]))
return true;
if (!std::isdigit(a[0]) && std::isdigit(b[0]))
return false;
if (!std::isdigit(a[0]) && !std::isdigit(b[0])) {
if (std::toupper(a[0]) == std::toupper(b[0]))
return compareNat(a.substr(1), b.substr(1));
return (std::toupper(a[0]) < std::toupper(b[0]));
}
// Both strings begin with digit --> parse both numbers
std::istringstream issa(a);
std::istringstream issb(b);
int ia, ib;
issa >> ia;
issb >> ib;
if (ia != ib)
return ia < ib;
// Numbers are the same --> remove numbers and recurse
std::string anew, bnew;
std::getline(issa, anew);
std::getline(issb, bnew);
return (compareNat(anew, bnew));
}
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;
}
inline std::string_view rightTrimSpaces(const std::string_view& str)
{
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;
}
inline std::string_view trimLeftAndRightSpaces(std::string_view str)
{
std::string_view trimmedString=str;
trimmedString = leftTrimSpaces(trimmedString);
trimmedString = rightTrimSpaces(trimmedString);
return trimmedString;
}
template<typename InputIt>
inline void normalize(InputIt itBegin, InputIt itEnd)
{
const auto squaredSumOfElements = std::accumulate(itBegin,
itEnd,
0.0,
[](const auto &sum, const auto &el) {
return sum + el * el;
});
assert(squaredSumOfElements != 0);
std::transform(itBegin, itEnd, itBegin, [&](auto &element) {
return element / std::sqrt(squaredSumOfElements);
});
}
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)
{
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;
@ -254,9 +149,8 @@ inline void parseIntegers(const std::string &str, std::vector<size_t> &result)
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]);
});
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) {
@ -285,6 +179,7 @@ inline std::vector<Vector6d> fromEigenMatrix(const Eigen::MatrixXd &m)
return v;
}
// std::string convertToLowercase(const std::string &s) {
// std::string lowercase;
// std::transform(s.begin(), s.end(), lowercase.begin(),
@ -311,19 +206,6 @@ inline std::vector<Vector6d> fromEigenMatrix(const Eigen::MatrixXd &m)
// return true;
//}
inline std::filesystem::path getFilepathWithExtension(const std::filesystem::path &folderPath,
const std::string &extension)
{
for (const std::filesystem::directory_entry &dirEntry :
std::filesystem::directory_iterator(folderPath)) {
if (dirEntry.is_regular_file() && std::filesystem::path(dirEntry).extension() == extension) {
return std::filesystem::path(dirEntry);
}
}
return "";
}
} // namespace Utilities
#ifdef POLYSCOPE_DEFINED
@ -344,7 +226,7 @@ inline void mainCallback()
// instead of full width. Must have
// matching PopItemWidth() below.
for (std::function<void()> &userCallback : globalPolyscopeData.userCallbacks) {
for (std::function<void()> userCallback : globalPolyscopeData.userCallbacks) {
userCallback();
}
@ -375,8 +257,7 @@ inline void init()
polyscope::view::upDir = polyscope::view::UpDir::ZUp;
polyscope::state::userCallback = &mainCallback;
polyscope::options::autocenterStructures = false;
polyscope::options::autoscaleStructures = false;
// polyscope::options::autocenterStructures = true;
}
using PolyscopeLabel = std::string;
inline std::pair<PolyscopeLabel, size_t> getSelection()
@ -394,12 +275,9 @@ inline void registerWorldAxes()
Eigen::MatrixX3d axesPositions(4, 3);
axesPositions.row(0) = Eigen::Vector3d(0, 0, 0);
// axesPositions.row(1) = Eigen::Vector3d(polyscope::state::lengthScale, 0, 0);
// axesPositions.row(2) = Eigen::Vector3d(0, polyscope::state::lengthScale, 0);
// axesPositions.row(3) = Eigen::Vector3d(0, 0, polyscope::state::lengthScale);
axesPositions.row(1) = Eigen::Vector3d(1, 0, 0);
axesPositions.row(2) = Eigen::Vector3d(0, 1, 0);
axesPositions.row(3) = Eigen::Vector3d(0, 0, 1);
axesPositions.row(1) = Eigen::Vector3d(polyscope::state::lengthScale, 0, 0);
axesPositions.row(2) = Eigen::Vector3d(0, polyscope::state::lengthScale, 0);
axesPositions.row(3) = Eigen::Vector3d(0, 0, polyscope::state::lengthScale);
Eigen::MatrixX2i axesEdges(3, 2);
axesEdges.row(0) = Eigen::Vector2i(0, 1);
@ -469,4 +347,17 @@ inline size_t computeHashOrdered(const std::vector<int> &v)
return std::hash<std::string>{}(elementsString);
}
//inline std::filesystem::path getFilepathWithExtension(const std::filesystem::path &folderPath,
// const std::string &extension)
//{
// for (const std::filesystem::directory_entry &dirEntry :
// std::filesystem::directory_iterator(folderPath)) {
// if (dirEntry.is_regular_file() && std::filesystem::path(dirEntry).extension() == extension) {
// return std::filesystem::path(dirEntry);
// }
// }
// return "";
//}
#endif // UTILITIES_H

8
vcgtrimesh.cpp
View File

@ -7,8 +7,8 @@
#include <wrap/io_trimesh/import.h>
//#include <wrap/nanoply/include/nanoplyWrapper.hpp>
bool VCGTriMesh::load(const std::filesystem::__cxx11::path &meshFilePath) {
assert(std::filesystem::exists(meshFilePath));
bool VCGTriMesh::load(const std::string &filename) {
assert(std::filesystem::exists(filename));
unsigned int mask = 0;
mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
@ -20,7 +20,7 @@ bool VCGTriMesh::load(const std::filesystem::__cxx11::path &meshFilePath) {
// std::filesystem::absolute(filename).string().c_str(), *this, mask)
// != 0) {
if (vcg::tri::io::Importer<VCGTriMesh>::Open(*this,
std::filesystem::absolute(meshFilePath).string().c_str())
std::filesystem::absolute(filename).string().c_str())
!= 0) {
std::cout << "Could not load tri mesh" << std::endl;
return false;
@ -32,7 +32,7 @@ bool VCGTriMesh::load(const std::filesystem::__cxx11::path &meshFilePath) {
vcg::tri::UpdateTopology<VCGTriMesh>::VertexEdge(*this);
vcg::tri::UpdateNormal<VCGTriMesh>::PerVertexNormalized(*this);
label = std::filesystem::path(meshFilePath).stem().string();
label = std::filesystem::path(filename).stem().string();
return true;
}

2
vcgtrimesh.hpp
View File

@ -39,7 +39,7 @@ class VCGTriMesh : public vcg::tri::TriMesh<std::vector<VCGTriMeshVertex>,
public:
VCGTriMesh();
VCGTriMesh(const std::string &filename);
bool load(const std::filesystem::path &meshFilePath) override;
bool load(const std::string &filename) override;
Eigen::MatrixX3d getVertices() const;
Eigen::MatrixX3i getFaces() const;
bool save(const std::string plyFilename);