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

Windows refactoring

This commit is contained in:
iasonmanolas 2021-12-23 21:11:27 +02:00
parent 31b8553b39 c9fc6ccd08
commit cd89112936
12 changed files with 642 additions and 485 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

89
CMakeLists.txt
View File

@ -1,4 +1,4 @@
cmake_minimum_required(VERSION 3.0) cmake_minimum_required(VERSION 2.8)
project(MySources) project(MySources)
set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_STANDARD_REQUIRED ON)
@ -9,30 +9,42 @@ if (CMAKE_VERSION VERSION_LESS 3.2)
else() else()
set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1") set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1")
endif() endif()
<<<<<<< HEAD
#set(EXTERNAL_DEPS_DIR "/home/iason/Coding/build/external dependencies") #set(EXTERNAL_DEPS_DIR "/home/iason/Coding/build/external dependencies")
if(NOT EXISTS EXTERNAL_DEPS_DIR) if(NOT EXISTS EXTERNAL_DEPS_DIR)
=======
set(EXTERNAL_DEPS_DIR "/home/iason/Coding/build/external dependencies")
if(NOT EXISTS ${EXTERNAL_DEPS_DIR})
>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb
set(EXTERNAL_DEPS_DIR ${CMAKE_CURRENT_BINARY_DIR}) set(EXTERNAL_DEPS_DIR ${CMAKE_CURRENT_BINARY_DIR})
message("External dependencies directory set:" ${EXTERNAL_DEPS_DIR}) message("External dependencies directory set:" ${EXTERNAL_DEPS_DIR})
endif() endif()
##Create directory for the external libraries ##Create directory for the external libraries
file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR}) file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR})
<<<<<<< HEAD
#message(${POLYSCOPE_ALREADY_COMPILED}) #message(${POLYSCOPE_ALREADY_COMPILED})
if(NOT DEFINED USE_POLYSCOPE) if(NOT DEFINED USE_POLYSCOPE)
message(FATAL_ERROR "Use polyscope was not defined") message(FATAL_ERROR "Use polyscope was not defined")
endif() endif()
=======
>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb
if(${USE_POLYSCOPE}) if(${USE_POLYSCOPE})
set(POLYSCOPE_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/polyscope)
download_project(PROJ POLYSCOPE download_project(PROJ POLYSCOPE
GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git
GIT_TAG master GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR} PREFIX ${EXTERNAL_DEPS_DIR}
BINARY_DIR ${POLYSCOPE_BINARY_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR}) add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR})
add_compile_definitions(POLYSCOPE_DEFINED)
endif() endif()
<<<<<<< HEAD
##dlib ##dlib
#set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin) #set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin)
@ -45,6 +57,8 @@ endif()
# ${UPDATE_DISCONNECTED_IF_AVAILABLE} # ${UPDATE_DISCONNECTED_IF_AVAILABLE}
#) #)
#add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR}) #add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR})
=======
>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb
##vcglib devel branch ##vcglib devel branch
download_project(PROJ vcglib_devel download_project(PROJ vcglib_devel
GIT_REPOSITORY https://github.com/IasonManolas/vcglib.git GIT_REPOSITORY https://github.com/IasonManolas/vcglib.git
@ -63,53 +77,80 @@ download_project(PROJ EIGEN
) )
add_subdirectory(${EIGEN_SOURCE_DIR} ${EIGEN_BINARY_DIR}) add_subdirectory(${EIGEN_SOURCE_DIR} ${EIGEN_BINARY_DIR})
##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 ##threed-beam-fea
set(threed-beam-fea_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/threed-beam-fea)
download_project(PROJ threed-beam-fea download_project(PROJ threed-beam-fea
GIT_REPOSITORY https://gitea-s2i2s.isti.cnr.it/manolas/threed-beam-fea.git GIT_REPOSITORY https://gitea-s2i2s.isti.cnr.it/manolas/threed-beam-fea.git
GIT_TAG master GIT_TAG master
BINARY_DIR ${threed-beam-fea_BINARY_DIR}
PREFIX ${EXTERNAL_DEPS_DIR} PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR}) add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR})
<<<<<<< HEAD
=======
##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)
>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb
if(MSVC) if(MSVC)
add_compile_definitions(_HAS_STD_BYTE=0) add_compile_definitions(_HAS_STD_BYTE=0)
endif(MSVC) endif(MSVC)
#link_directories(${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) 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) add_library(${PROJECT_NAME} ${MySourcesFiles} ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)
target_include_directories(${PROJECT_NAME}
set(USE_TBB true) PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph
if(${USE_TBB}) PUBLIC ${vcglib_devel_SOURCE_DIR}
##TBB PUBLIC ${threed-beam-fea_SOURCE_DIR}
download_project(PROJ TBB )
GIT_REPOSITORY https://github.com/wjakob/tbb.git if(${MYSOURCES_STATIC_LINK})
GIT_TAG master message("Linking statically")
PREFIX ${EXTERNAL_DEPS_DIR} target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread gfortran quadmath)
${UPDATE_DISCONNECTED_IF_AVAILABLE} else()
) target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot ThreedBeamFEA tbb pthread)
option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF) if(${USE_POLYSCOPE})
option(TBB_BUILD_STATIC "Build TBB static library" ON) target_link_libraries(${PROJECT_NAME} polyscope)
add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR})
link_directories(${TBB_BINARY_DIR})
if(${MYSOURCES_STATIC_LINK})
target_link_libraries(${PROJECT_NAME} PUBLIC -static tbb_static)
else()
target_link_libraries(${PROJECT_NAME} PUBLIC tbb)
endif() endif()
endif() endif()
target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
if(USE_ENSMALLEN) if(USE_ENSMALLEN)
##ENSMALLEN ##ENSMALLEN
set(ENSMALLEN_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/ensmallen)
download_project(PROJ ENSMALLEN download_project(PROJ ENSMALLEN
GIT_REPOSITORY https://github.com/mlpack/ensmallen.git GIT_REPOSITORY https://github.com/mlpack/ensmallen.git
GIT_TAG master GIT_TAG master
BINARY_DIR ${ENSMALLEN_BINARY_DIR}
PREFIX ${EXTERNAL_DEPS_DIR} PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
add_subdirectory(${ENSMALLEN_SOURCE_DIR} ${ENSMALLEN_BINARY_DIR}) add_subdirectory(${ENSMALLEN_SOURCE_DIR} ${ENSMALLEN_BINARY_DIR})
<<<<<<< HEAD
download_project(PROJ ARMADILLO download_project(PROJ ARMADILLO
GIT_REPOSITORY https://gitlab.com/conradsnicta/armadillo-code.git GIT_REPOSITORY https://gitlab.com/conradsnicta/armadillo-code.git
@ -120,14 +161,18 @@ if(USE_ENSMALLEN)
add_subdirectory(${ARMADILLO_SOURCE_DIR} ${ARMADILLO_BINARY_DIR}) add_subdirectory(${ARMADILLO_SOURCE_DIR} ${ARMADILLO_BINARY_DIR})
target_include_directories(${PROJECT_NAME} PUBLIC armadillo ) target_include_directories(${PROJECT_NAME} PUBLIC armadillo )
=======
set(ARMADILLO_SOURCE_DIR "/home/iason/Coding/Libraries/armadillo")
add_subdirectory(${ARMADILLO_SOURCE_DIR} ${EXTERNAL_DEPS_DIR}/armadillo_build)
>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb
if(${MYSOURCES_STATIC_LINK}) if(${MYSOURCES_STATIC_LINK})
target_link_libraries(${PROJECT_NAME} PUBLIC "-static" ensmallen ${EXTERNAL_DEPS_DIR}/armadillo_build/libarmadillo.a) target_link_libraries(${PROJECT_NAME} PUBLIC "-static" ensmallen ${EXTERNAL_DEPS_DIR}/armadillo_build/libarmadillo.a)
else() else()
target_link_libraries(${PROJECT_NAME} PUBLIC armadillo ensmallen) target_link_libraries(${PROJECT_NAME} PUBLIC armadillo ensmallen)
endif() endif()
target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include) target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include)
endif() endif()
<<<<<<< HEAD
target_include_directories(${PROJECT_NAME} target_include_directories(${PROJECT_NAME}
PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph
@ -183,3 +228,5 @@ endif()
#target_link_libraries(${PROJECT_NAME} PUBLIC GSL) #target_link_libraries(${PROJECT_NAME} PUBLIC GSL)
#target_include_directories(${PROJECT_NAME} PUBLIC GSL) #target_include_directories(${PROJECT_NAME} PUBLIC GSL)
=======
>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb

34
beam.hpp
View File

@ -9,15 +9,36 @@ struct RectangularBeamDimensions {
inline static std::string name{"Rectangular"}; inline static std::string name{"Rectangular"};
double b; double b;
double h; double h;
RectangularBeamDimensions(const double &width, const double &height)
: b(width), h(height) { double A{0}; // cross sectional area
assert(width > 0 && height > 0);
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)
{
assert(width > 0 && height > 0);
updateProperties();
} }
RectangularBeamDimensions() : b(0.002), h(0.002) {} RectangularBeamDimensions() : b(0.002), h(0.002) { updateProperties(); }
std::string toString() const { std::string toString() const {
return std::string("b=") + std::to_string(b) + std::string(" h=") + return std::string("b=") + std::to_string(b) + std::string(" h=") +
std::to_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 { struct CylindricalBeamDimensions {
@ -39,18 +60,21 @@ struct ElementMaterial
{ {
double poissonsRatio; double poissonsRatio;
double youngsModulus; double youngsModulus;
double G;
ElementMaterial(const double &poissonsRatio, const double &youngsModulus) ElementMaterial(const double &poissonsRatio, const double &youngsModulus)
: poissonsRatio(poissonsRatio), youngsModulus(youngsModulus) : poissonsRatio(poissonsRatio), youngsModulus(youngsModulus)
{ {
assert(poissonsRatio <= 0.5 && poissonsRatio >= -1); assert(poissonsRatio <= 0.5 && poissonsRatio >= -1);
updateProperties();
} }
ElementMaterial() : poissonsRatio(0.3), youngsModulus(200 * 1e9) {} ElementMaterial() : poissonsRatio(0.3), youngsModulus(200 * 1e9) { updateProperties(); }
std::string toString() const std::string toString() const
{ {
return std::string("Material:") + std::string("\nPoisson's ratio=") return std::string("Material:") + std::string("\nPoisson's ratio=")
+ std::to_string(poissonsRatio) + std::string("\nYoung's Modulus(GPa)=") + std::to_string(poissonsRatio) + std::string("\nYoung's Modulus(GPa)=")
+ std::to_string(youngsModulus / 1e9); + std::to_string(youngsModulus / 1e9);
} }
void updateProperties() { G = youngsModulus / (2 * (1 + poissonsRatio)); }
}; };
#endif // BEAM_HPP #endif // BEAM_HPP

41
drmsimulationmodel.cpp
View File

@ -796,21 +796,21 @@ double DRMSimulationModel::computeTotalInternalPotentialEnergy()
const EdgeIndex ei = pMesh->getIndex(e); const EdgeIndex ei = pMesh->getIndex(e);
const double e_k = element.length - element.initialLength; const double e_k = element.length - element.initialLength;
const double axialPE = pow(e_k, 2) * element.material.youngsModulus * element.A const double axialPE = pow(e_k, 2) * element.material.youngsModulus * element.dimensions.A
/ (2 * element.initialLength); / (2 * element.initialLength);
const double theta1Diff = theta1_jplus1 - theta1_j; const double theta1Diff = theta1_jplus1 - theta1_j;
const double torsionalPE = element.G * element.inertia.J * pow(theta1Diff, 2) const double torsionalPE = element.material.G * element.dimensions.inertia.J
/ (2 * element.initialLength); * pow(theta1Diff, 2) / (2 * element.initialLength);
const double firstBendingPE_inBracketsTerm = 4 * pow(theta2_j, 2) const double firstBendingPE_inBracketsTerm = 4 * pow(theta2_j, 2)
+ 4 * theta2_j * theta2_jplus1 + 4 * theta2_j * theta2_jplus1
+ 4 * pow(theta2_jplus1, 2); + 4 * pow(theta2_jplus1, 2);
const double firstBendingPE = firstBendingPE_inBracketsTerm * element.material.youngsModulus const double firstBendingPE = firstBendingPE_inBracketsTerm * element.material.youngsModulus
* element.inertia.I2 / (2 * element.initialLength); * element.dimensions.inertia.I2 / (2 * element.initialLength);
const double secondBendingPE_inBracketsTerm = 4 * pow(theta3_j, 2) const double secondBendingPE_inBracketsTerm = 4 * pow(theta3_j, 2)
+ 4 * theta3_j * theta3_jplus1 + 4 * theta3_j * theta3_jplus1
+ 4 * pow(theta3_jplus1, 2); + 4 * pow(theta3_jplus1, 2);
const double secondBendingPE = secondBendingPE_inBracketsTerm const double secondBendingPE = secondBendingPE_inBracketsTerm
* element.material.youngsModulus * element.inertia.I3 * element.material.youngsModulus * element.dimensions.inertia.I3
/ (2 * element.initialLength); / (2 * element.initialLength);
totalInternalPotentialEnergy += axialPE + torsionalPE + firstBendingPE + secondBendingPE; totalInternalPotentialEnergy += axialPE + torsionalPE + firstBendingPE + secondBendingPE;
@ -1261,14 +1261,15 @@ void DRMSimulationModel::updateNodalMasses()
for (const EdgePointer &ep : pMesh->nodes[v].incidentElements) { for (const EdgePointer &ep : pMesh->nodes[v].incidentElements) {
const size_t ei = pMesh->getIndex(ep); const size_t ei = pMesh->getIndex(ep);
const Element &elem = pMesh->elements[ep]; const Element &elem = pMesh->elements[ep];
const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length; const double SkTranslational = elem.material.youngsModulus * elem.dimensions.A
/ elem.length;
translationalSumSk += SkTranslational; translationalSumSk += SkTranslational;
const double lengthToThe3 = std::pow(elem.length, 3); const double lengthToThe3 = std::pow(elem.length, 3);
const long double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2 const long double SkRotational_I2 = elem.material.youngsModulus * elem.dimensions.inertia.I2
/ lengthToThe3; / lengthToThe3;
const long double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3 const long double SkRotational_I3 = elem.material.youngsModulus * elem.dimensions.inertia.I3
/ lengthToThe3; / lengthToThe3;
const long double SkRotational_J = elem.material.youngsModulus * elem.inertia.J const long double SkRotational_J = elem.material.youngsModulus * elem.dimensions.inertia.J
/ lengthToThe3; / lengthToThe3;
rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I2 += SkRotational_I2;
rotationalSumSk_I3 += SkRotational_I3; rotationalSumSk_I3 += SkRotational_I3;
@ -1627,14 +1628,15 @@ void DRMSimulationModel::updatePositionsOnTheFly(
double rotationalSumSk_J = 0; double rotationalSumSk_J = 0;
for (const EdgePointer &ep : pMesh->nodes[v].incidentElements) { for (const EdgePointer &ep : pMesh->nodes[v].incidentElements) {
const Element &elem = pMesh->elements[ep]; const Element &elem = pMesh->elements[ep];
const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length; const double SkTranslational = elem.material.youngsModulus * elem.dimensions.A
/ elem.length;
translationalSumSk += SkTranslational; translationalSumSk += SkTranslational;
const double lengthToThe3 = std::pow(elem.length, 3); const double lengthToThe3 = std::pow(elem.length, 3);
const double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2 const double SkRotational_I2 = elem.material.youngsModulus * elem.dimensions.inertia.I2
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3 const double SkRotational_I3 = elem.material.youngsModulus * elem.dimensions.inertia.I3
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const double SkRotational_J = elem.material.youngsModulus * elem.inertia.J const double SkRotational_J = elem.material.youngsModulus * elem.dimensions.inertia.J
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I2 += SkRotational_I2;
rotationalSumSk_I3 += SkRotational_I3; rotationalSumSk_I3 += SkRotational_I3;
@ -1831,10 +1833,13 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptr<Simul
results.debug_q_normal[vi] = q_normal; results.debug_q_normal[vi] = q_normal;
results.debug_q_nr[vi] = q_nr_nInit; results.debug_q_nr[vi] = q_nr_nInit;
results.rotationalDisplacementQuaternion[vi] results.rotationalDisplacementQuaternion[vi]
//Eigen::Quaterniond R
= (q_normal = (q_normal
* (q_f1_nInit * q_nr_nInit)); //q_f1_nDeformed * q_nr_nDeformed * q_normal; * (q_f1_nInit * q_nr_nInit)); //q_f1_nDeformed * q_nr_nDeformed * q_normal;
//Update the displacement vector to contain the euler angles //Update the displacement vector to contain the euler angles
const Eigen::Vector3d eulerAngles = results.rotationalDisplacementQuaternion[vi] const Eigen::Vector3d eulerAngles = results
.rotationalDisplacementQuaternion[vi]
// R
.toRotationMatrix() .toRotationMatrix()
.eulerAngles(0, 1, 2); .eulerAngles(0, 1, 2);
results.displacements[vi][3] = eulerAngles[0]; results.displacements[vi][3] = eulerAngles[0];
@ -2007,10 +2012,10 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
std::vector<double> I3(pMesh->EN()); std::vector<double> I3(pMesh->EN());
for (const EdgeType &e : pMesh->edge) { for (const EdgeType &e : pMesh->edge) {
const size_t ei = pMesh->getIndex(e); const size_t ei = pMesh->getIndex(e);
A[ei] = pMesh->elements[e].A; A[ei] = pMesh->elements[e].dimensions.A;
J[ei] = pMesh->elements[e].inertia.J; J[ei] = pMesh->elements[e].dimensions.inertia.J;
I2[ei] = pMesh->elements[e].inertia.I2; I2[ei] = pMesh->elements[e].dimensions.inertia.I2;
I3[ei] = pMesh->elements[e].inertia.I3; I3[ei] = pMesh->elements[e].dimensions.inertia.I3;
} }
polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A); polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A);

3
drmsimulationmodel.hpp
View File

@ -2,7 +2,8 @@
#define BEAMFORMFINDER_HPP #define BEAMFORMFINDER_HPP
//#ifdef USE_MATPLOT //#ifdef USE_MATPLOT
#include "matplot/matplot.h" //#include "matplot.h"
#include <matplot/matplot.h>
//#endif //#endif
#include "simulationmesh.hpp" #include "simulationmesh.hpp"
#include "simulationmodel.hpp" #include "simulationmodel.hpp"

40
linearsimulationmodel.cpp
View File

@ -1,5 +1,5 @@
#include "linearsimulationmodel.hpp" #include "linearsimulationmodel.hpp"
#include "gsl/gsl" //#include "gsl/gsl"
std::vector<fea::Elem> LinearSimulationModel::getFeaElements( std::vector<fea::Elem> LinearSimulationModel::getFeaElements(
const std::shared_ptr<SimulationMesh> &pMesh) const std::shared_ptr<SimulationMesh> &pMesh)
@ -14,10 +14,10 @@ std::vector<fea::Elem> LinearSimulationModel::getFeaElements(
(evn0[2] + evn1[2]) / 2}; (evn0[2] + evn1[2]) / 2};
const Element &element = pMesh->elements[edgeIndex]; const Element &element = pMesh->elements[edgeIndex];
const double E = element.material.youngsModulus; const double E = element.material.youngsModulus;
fea::Props feaProperties(E * element.A, fea::Props feaProperties(E * element.dimensions.A,
E * element.inertia.I3, E * element.dimensions.inertia.I3,
E * element.inertia.I2, E * element.dimensions.inertia.I2,
element.G * element.inertia.J, element.material.G * element.dimensions.inertia.J,
nAverageVector); nAverageVector);
const int vi0 = pMesh->getIndex(pMesh->edge[edgeIndex].cV(0)); const int vi0 = pMesh->getIndex(pMesh->edge[edgeIndex].cV(0));
const int vi1 = pMesh->getIndex(pMesh->edge[edgeIndex].cV(1)); const int vi1 = pMesh->getIndex(pMesh->edge[edgeIndex].cV(1));
@ -139,26 +139,30 @@ SimulationResults LinearSimulationModel::getResults(
const double elementPotentialEnergy_axial_v0 = std::pow(elementForces[0], 2) const double elementPotentialEnergy_axial_v0 = std::pow(elementForces[0], 2)
* element.initialLength * element.initialLength
/ (4 * element.material.youngsModulus / (4 * element.material.youngsModulus
* element.A); * element.dimensions.A);
const double elementPotentialEnergy_axial_v1 = std::pow(elementForces[6], 2) const double elementPotentialEnergy_axial_v1 = std::pow(elementForces[6], 2)
* element.initialLength * element.initialLength
/ (4 * element.material.youngsModulus / (4 * element.material.youngsModulus
* element.A); * element.dimensions.A);
const double elementPotentialEnergy_axial = elementPotentialEnergy_axial_v0 const double elementPotentialEnergy_axial = elementPotentialEnergy_axial_v0
+ elementPotentialEnergy_axial_v1; + elementPotentialEnergy_axial_v1;
//Shear //Shear
const double elementPotentialEnergy_shearY_v0 = std::pow(elementForces[1], 2) const double elementPotentialEnergy_shearY_v0 = std::pow(elementForces[1], 2)
* element.initialLength * element.initialLength
/ (4 * element.A * element.G); / (4 * element.dimensions.A
* element.material.G);
const double elementPotentialEnergy_shearZ_v0 = std::pow(elementForces[2], 2) const double elementPotentialEnergy_shearZ_v0 = std::pow(elementForces[2], 2)
* element.initialLength * element.initialLength
/ (4 * element.A * element.G); / (4 * element.dimensions.A
* element.material.G);
const double elementPotentialEnergy_shearY_v1 = std::pow(elementForces[7], 2) const double elementPotentialEnergy_shearY_v1 = std::pow(elementForces[7], 2)
* element.initialLength * element.initialLength
/ (4 * element.A * element.G); / (4 * element.dimensions.A
* element.material.G);
const double elementPotentialEnergy_shearZ_v1 = std::pow(elementForces[8], 2) const double elementPotentialEnergy_shearZ_v1 = std::pow(elementForces[8], 2)
* element.initialLength * element.initialLength
/ (4 * element.A * element.G); / (4 * element.dimensions.A
* element.material.G);
const double elementPotentialEnergy_shear = elementPotentialEnergy_shearY_v0 const double elementPotentialEnergy_shear = elementPotentialEnergy_shearY_v0
+ elementPotentialEnergy_shearZ_v0 + elementPotentialEnergy_shearZ_v0
+ elementPotentialEnergy_shearY_v1 + elementPotentialEnergy_shearY_v1
@ -167,19 +171,19 @@ SimulationResults LinearSimulationModel::getResults(
const double elementPotentialEnergy_bendingY_v0 = std::pow(elementForces[4], 2) const double elementPotentialEnergy_bendingY_v0 = std::pow(elementForces[4], 2)
* element.initialLength * element.initialLength
/ (4 * element.material.youngsModulus / (4 * element.material.youngsModulus
* element.inertia.I2); * element.dimensions.inertia.I2);
const double elementPotentialEnergy_bendingZ_v0 = std::pow(elementForces[5], 2) const double elementPotentialEnergy_bendingZ_v0 = std::pow(elementForces[5], 2)
* element.initialLength * element.initialLength
/ (4 * element.material.youngsModulus / (4 * element.material.youngsModulus
* element.inertia.I3); * element.dimensions.inertia.I3);
const double elementPotentialEnergy_bendingY_v1 = std::pow(elementForces[10], 2) const double elementPotentialEnergy_bendingY_v1 = std::pow(elementForces[10], 2)
* element.initialLength * element.initialLength
/ (4 * element.material.youngsModulus / (4 * element.material.youngsModulus
* element.inertia.I2); * element.dimensions.inertia.I2);
const double elementPotentialEnergy_bendingZ_v1 = std::pow(elementForces[11], 2) const double elementPotentialEnergy_bendingZ_v1 = std::pow(elementForces[11], 2)
* element.initialLength * element.initialLength
/ (4 * element.material.youngsModulus / (4 * element.material.youngsModulus
* element.inertia.I3); * element.dimensions.inertia.I3);
const double elementPotentialEnergy_bending = elementPotentialEnergy_bendingY_v0 const double elementPotentialEnergy_bending = elementPotentialEnergy_bendingY_v0
+ elementPotentialEnergy_bendingZ_v0 + elementPotentialEnergy_bendingZ_v0
+ elementPotentialEnergy_bendingY_v1 + elementPotentialEnergy_bendingY_v1
@ -187,10 +191,12 @@ SimulationResults LinearSimulationModel::getResults(
//Torsion //Torsion
const double elementPotentialEnergy_torsion_v0 = std::pow(elementForces[3], 2) const double elementPotentialEnergy_torsion_v0 = std::pow(elementForces[3], 2)
* element.initialLength * element.initialLength
/ (4 * element.inertia.J * element.G); / (4 * element.dimensions.inertia.J
* element.material.G);
const double elementPotentialEnergy_torsion_v1 = std::pow(elementForces[9], 2) const double elementPotentialEnergy_torsion_v1 = std::pow(elementForces[9], 2)
* element.initialLength * element.initialLength
/ (4 * element.inertia.J * element.G); / (4 * element.dimensions.inertia.J
* element.material.G);
const double elementPotentialEnergy_torsion = elementPotentialEnergy_torsion_v0 const double elementPotentialEnergy_torsion = elementPotentialEnergy_torsion_v0
+ elementPotentialEnergy_torsion_v1; + elementPotentialEnergy_torsion_v1;
const double elementInternalPotentialEnergy = elementPotentialEnergy_axial const double elementInternalPotentialEnergy = elementPotentialEnergy_axial

2
mesh.hpp
View File

@ -1,7 +1,7 @@
#ifndef MESH_HPP #ifndef MESH_HPP
#define MESH_HPP #define MESH_HPP
#include <Eigen/Core> //#include <Eigen/Core>
#include <string> #include <string>
class Mesh class Mesh

275
reducedmodeloptimizer_structs.hpp
View File

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

319
simulation_structs.hpp
View File

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

134
simulationhistoryplotter.hpp
View File

@ -114,93 +114,81 @@ struct SimulationResultsReporter {
if (YvaluesToPlot.size() < 2) { if (YvaluesToPlot.size() < 2) {
return; return;
} }
auto x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size()); std::vector<double> colors(YvaluesToPlot.size(), 0.5);
std::vector<double> colors(x.size(), 0.5); std::vector<double> markerSizes(YvaluesToPlot.size(), 5);
std::vector<double> markerSizes(x.size(), 5);
if (!markPoints.empty()) { if (!markPoints.empty()) {
for (const auto pointIndex : markPoints) { for (const auto pointIndex : markPoints) {
colors[pointIndex] = 0.9; colors[pointIndex] = 0.9;
markerSizes[pointIndex] = 14; markerSizes[pointIndex] = 14;
} }
} }
std::vector<double> x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size());
createPlot(xLabel, yLabel, x, YvaluesToPlot, markerSizes, colors, saveTo); 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, void createPlots(const SimulationHistory &history,
const std::string &reportFolderPath, const std::string &reportFolderPath,
const std::string &graphSuffix) { const std::string &graphSuffix)
const auto graphsFolder = {
std::filesystem::path(reportFolderPath).append("Graphs"); const auto graphsFolder = std::filesystem::path(reportFolderPath).append("Graphs");
std::filesystem::remove_all(graphsFolder); std::filesystem::remove_all(graphsFolder);
std::filesystem::create_directory(graphsFolder.string()); std::filesystem::create_directory(graphsFolder.string());
if (!history.kineticEnergy.empty()) { if (!history.kineticEnergy.empty()) {
createPlot("Number of Iterations", createPlot("Number of Iterations",
"Log of Kinetic Energy log", "Log of Kinetic Energy log",
history.kineticEnergy, history.kineticEnergy,
std::filesystem::path(graphsFolder) std::filesystem::path(graphsFolder)
.append("KineticEnergyLog_" + graphSuffix + ".png") .append("KineticEnergyLog_" + graphSuffix + ".png")
.string(), .string(),
history.redMarks); history.redMarks);
} }
if (!history.logResidualForces.empty()) { if (!history.logResidualForces.empty()) {
createPlot("Number of Iterations", createPlot("Number of Iterations",
"Residual Forces norm log", "Residual Forces norm log",
history.logResidualForces, history.logResidualForces,
std::filesystem::path(graphsFolder) std::filesystem::path(graphsFolder)
.append("ResidualForcesLog_" + graphSuffix + ".png") .append("ResidualForcesLog_" + graphSuffix + ".png")
.string(), .string(),
history.redMarks); history.redMarks);
} }
if (!history.potentialEnergies.empty()) { if (!history.potentialEnergies.empty()) {
createPlot("Number of Iterations", createPlot("Number of Iterations",
"Potential energy", "Potential energy",
history.potentialEnergies, history.potentialEnergies,
std::filesystem::path(graphsFolder) std::filesystem::path(graphsFolder)
.append("PotentialEnergy_" + graphSuffix + ".png") .append("PotentialEnergy_" + graphSuffix + ".png")
.string(), .string(),
history.redMarks); history.redMarks);
} }
if (!history.residualForcesMovingAverage.empty()) { if (!history.residualForcesMovingAverage.empty()) {
createPlot("Number of Iterations", createPlot("Number of Iterations",
"Residual forces moving average", "Residual forces moving average",
history.residualForcesMovingAverage, history.residualForcesMovingAverage,
std::filesystem::path(graphsFolder) std::filesystem::path(graphsFolder)
.append("ResidualForcesMovingAverage_" + graphSuffix + ".png") .append("ResidualForcesMovingAverage_" + graphSuffix + ".png")
.string(), .string(),
history.redMarks); history.redMarks);
} }
// if (!history.residualForcesMovingAverageDerivativesLog.empty()) { // if (!history.residualForcesMovingAverageDerivativesLog.empty()) {
// createPlot("Number of Iterations", // createPlot("Number of Iterations",
// "Residual forces moving average derivative log", // "Residual forces moving average derivative log",
// history.residualForcesMovingAverageDerivativesLog, // history.residualForcesMovingAverageDerivativesLog,
// std::filesystem::path(graphsFolder) // std::filesystem::path(graphsFolder)
// .append("ResidualForcesMovingAverageDerivativesLog_" + graphSuffix + ".png") // .append("ResidualForcesMovingAverageDerivativesLog_" + graphSuffix + ".png")
// .string()); // .string());
// } // }
if (!history.sumOfNormalizedDisplacementNorms.empty()) { if (!history.sumOfNormalizedDisplacementNorms.empty()) {
createPlot("Number of Iterations", createPlot("Number of Iterations",
"Sum of normalized displacement norms", "Sum of normalized displacement norms",
history.sumOfNormalizedDisplacementNorms, history.sumOfNormalizedDisplacementNorms,
std::filesystem::path(graphsFolder) std::filesystem::path(graphsFolder)
.append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png") .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png")
.string(), .string(),
history.redMarks); history.redMarks);
} }
} }
}; };

79
simulationmesh.cpp
View File

@ -198,7 +198,6 @@ void SimulationMesh::reset() {
} }
void SimulationMesh::initializeElements() { void SimulationMesh::initializeElements() {
computeElementalProperties();
for (const EdgeType &e : edge) { for (const EdgeType &e : edge) {
Element &element = elements[e]; Element &element = elements[e];
element.ei = getIndex(e); element.ei = getIndex(e);
@ -258,7 +257,6 @@ void SimulationMesh::setBeamCrossSection(
const CrossSectionType &beamDimensions) { const CrossSectionType &beamDimensions) {
for (size_t ei = 0; ei < EN(); ei++) { for (size_t ei = 0; ei < EN(); ei++) {
elements[ei].dimensions = beamDimensions; elements[ei].dimensions = beamDimensions;
elements[ei].computeDimensionsProperties(beamDimensions);
elements[ei].updateRigidity(); elements[ei].updateRigidity();
} }
} }
@ -455,60 +453,59 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1,
return angle; return angle;
} }
void Element::computeMaterialProperties(const ElementMaterial &material) { //void Element::computeMaterialProperties(const ElementMaterial &material) {
G = material.youngsModulus / (2 * (1 + material.poissonsRatio)); // G = material.youngsModulus / (2 * (1 + material.poissonsRatio));
} //}
void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A) //void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A)
{ //{
A = dimensions.b * dimensions.h; // A = dimensions.b * dimensions.h;
} //}
void Element::computeDimensionsProperties( //void Element::computeDimensionsProperties(
const RectangularBeamDimensions &dimensions) { // const RectangularBeamDimensions &dimensions) {
assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions)); // assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions));
computeCrossSectionArea(dimensions, A); // computeCrossSectionArea(dimensions, A);
computeMomentsOfInertia(dimensions, inertia); // computeMomentsOfInertia(dimensions, dimensions.inertia);
} //}
void Element::computeDimensionsProperties( //void Element::computeDimensionsProperties(
const CylindricalBeamDimensions &dimensions) { // const CylindricalBeamDimensions &dimensions) {
assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); // assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2)); // 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; // dimensions.inertia.I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64;
inertia.I3 = inertia.I2; // dimensions.inertia.I3 = dimensions.inertia.I2;
inertia.J = inertia.I2 + inertia.I3; // dimensions.inertia.J = dimensions.inertia.I2 + dimensions.inertia.I3;
} //}
void Element::setDimensions(const CrossSectionType &dimensions) { void Element::setDimensions(const CrossSectionType &dimensions) {
this->dimensions = dimensions; this->dimensions = dimensions;
computeDimensionsProperties(dimensions); assert(this->dimensions.A == dimensions.A);
// computeDimensionsProperties(dimensions);
updateRigidity(); updateRigidity();
} }
void Element::setMaterial(const ElementMaterial &material) { void Element::setMaterial(const ElementMaterial &material)
this->material = material; {
computeMaterialProperties(material); this->material = material;
updateRigidity(); // computeMaterialProperties(material);
updateRigidity();
} }
double Element::getMass(const double &materialDensity) double Element::getMass(const double &materialDensity)
{ {
const double beamVolume = A * length; const double beamVolume = dimensions.A * length;
return beamVolume * materialDensity; return beamVolume * materialDensity;
} }
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() { void Element::updateRigidity() {
rigidity.axial = material.youngsModulus * A / initialLength; // assert(initialLength != 0);
rigidity.torsional = G * inertia.J / initialLength; rigidity.axial = material.youngsModulus * dimensions.A / initialLength;
rigidity.firstBending = 2 * material.youngsModulus * inertia.I2 / initialLength; // assert(rigidity.axial != 0);
rigidity.secondBending = 2 * material.youngsModulus * inertia.I3 / initialLength; 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);
} }

18
simulationmesh.hpp
View File

@ -1,9 +1,9 @@
#ifndef SIMULATIONMESH_HPP #ifndef SIMULATIONMESH_HPP
#define SIMULATIONMESH_HPP #define SIMULATIONMESH_HPP
#include "Eigen/Dense"
#include "edgemesh.hpp" #include "edgemesh.hpp"
#include "trianglepatterngeometry.hpp" #include "trianglepatterngeometry.hpp"
#include <Eigen/Dense>
//extern bool drawGlobal; //extern bool drawGlobal;
struct Element; struct Element;
@ -62,26 +62,14 @@ struct Element {
CrossSectionType dimensions; CrossSectionType dimensions;
ElementMaterial material; ElementMaterial material;
double G{0};
double A{0}; // cross sectional area
void computeMaterialProperties(const ElementMaterial &material); void computeMaterialProperties(const ElementMaterial &material);
void computeDimensionsProperties(const RectangularBeamDimensions &dimensions); // void computeDimensionsProperties(const RectangularBeamDimensions &dimensions);
void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions); // void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions);
void setDimensions(const CrossSectionType &dimensions); void setDimensions(const CrossSectionType &dimensions);
void setMaterial(const ElementMaterial &material); void setMaterial(const ElementMaterial &material);
double getMass(const double &matrialDensity); 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 { struct LocalFrame {
VectorType t1; VectorType t1;
VectorType t2; VectorType t2;

93
utilities.hpp
View File

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