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

Working refactoring

This commit is contained in:
iasonmanolas 2021-12-19 20:15:36 +02:00
parent 1966207b4c
commit 56e5e043f6
11 changed files with 238 additions and 284 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

156
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,7 +9,6 @@ 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()
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 ${CMAKE_CURRENT_BINARY_DIR}) set(EXTERNAL_DEPS_DIR ${CMAKE_CURRENT_BINARY_DIR})
@ -17,30 +16,27 @@ endif()
##Create directory for the external libraries ##Create directory for the external libraries
file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR}) file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR})
#message(${POLYSCOPE_ALREADY_COMPILED}) #message(${POLYSCOPE_ALREADY_COMPILED})
if(${USE_POLYSCOPE}) if(${USE_POLYSCOPE})
download_project(PROJ POLYSCOPE download_project(PROJ POLYSCOPE
GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git
GIT_TAG master GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR} PREFIX ${EXTERNAL_DEPS_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()
#dlib
set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin)
##dlib file(MAKE_DIRECTORY ${DLIB_BIN_DIR})
#set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin) download_project(PROJ DLIB
#file(MAKE_DIRECTORY ${DLIB_BIN_DIR}) GIT_REPOSITORY https://github.com/davisking/dlib.git
#download_project(PROJ DLIB GIT_TAG master
# GIT_REPOSITORY https://github.com/davisking/dlib.git BINARY_DIR ${DLIB_BIN_DIR}
# GIT_TAG master PREFIX ${EXTERNAL_DEPS_DIR}
# BINARY_DIR ${DLIB_BIN_DIR} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
# PREFIX ${EXTERNAL_DEPS_DIR} )
# ${UPDATE_DISCONNECTED_IF_AVAILABLE} add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR})
#)
#add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR})
##vcglib devel branch ##vcglib devel branch
download_project(PROJ vcglib_devel download_project(PROJ vcglib_devel
GIT_REPOSITORY https://github.com/IasonManolas/vcglib.git GIT_REPOSITORY https://github.com/IasonManolas/vcglib.git
@ -49,7 +45,14 @@ download_project(PROJ vcglib_devel
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
add_subdirectory(${vcglib_devel_SOURCE_DIR} ${vcglib_devel_BINARY_DIR}) add_subdirectory(${vcglib_devel_SOURCE_DIR} ${vcglib_devel_BINARY_DIR})
##matplot++ lib
download_project(PROJ MATPLOTPLUSPLUS
GIT_REPOSITORY https://github.com/alandefreitas/matplotplusplus
GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${MATPLOTPLUSPLUS_SOURCE_DIR} ${MATPLOTPLUSPLUS_BINARY_DIR})
##threed-beam-fea ##threed-beam-fea
download_project(PROJ threed-beam-fea download_project(PROJ threed-beam-fea
GIT_REPOSITORY https://github.com/IasonManolas/threed-beam-fea.git GIT_REPOSITORY https://github.com/IasonManolas/threed-beam-fea.git
@ -58,37 +61,40 @@ download_project(PROJ threed-beam-fea
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR}) add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR})
##TBB
download_project(PROJ TBB
GIT_REPOSITORY https://github.com/wjakob/tbb.git
GIT_TAG master
PREFIX ${EXTERNAL_DEPS_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 ###Eigen 3 NOTE: Eigen is required on the system the code is ran
find_package(Eigen3 3.3 REQUIRED) find_package(Eigen3 3.3 REQUIRED)
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 # PUBLIC ${MySourcesFiles}
GIT_REPOSITORY https://github.com/wjakob/tbb.git )
GIT_TAG master if(${MYSOURCES_STATIC_LINK})
PREFIX ${EXTERNAL_DEPS_DIR} message("Linking statically")
${UPDATE_DISCONNECTED_IF_AVAILABLE} target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread gfortran quadmath)
) else()
option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF) target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib ThreedBeamFEA tbb pthread)
option(TBB_BUILD_STATIC "Build TBB static library" ON) if(${USE_POLYSCOPE})
add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR}) target_link_libraries(${PROJECT_NAME} polyscope)
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
@ -99,8 +105,6 @@ if(USE_ENSMALLEN)
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
add_subdirectory(${ENSMALLEN_SOURCE_DIR} ${ENSMALLEN_BINARY_DIR}) add_subdirectory(${ENSMALLEN_SOURCE_DIR} ${ENSMALLEN_BINARY_DIR})
set(ARMADILLO_SOURCE_DIR "/home/iason/Coding/Libraries/armadillo") set(ARMADILLO_SOURCE_DIR "/home/iason/Coding/Libraries/armadillo")
add_subdirectory(${ARMADILLO_SOURCE_DIR} ${EXTERNAL_DEPS_DIR}/armadillo_build) add_subdirectory(${ARMADILLO_SOURCE_DIR} ${EXTERNAL_DEPS_DIR}/armadillo_build)
if(${MYSOURCES_STATIC_LINK}) if(${MYSOURCES_STATIC_LINK})
@ -109,61 +113,5 @@ if(USE_ENSMALLEN)
target_link_libraries(${PROJECT_NAME} PUBLIC) target_link_libraries(${PROJECT_NAME} PUBLIC)
target_link_libraries(${PROJECT_NAME} PUBLIC armadillo ensmallen) target_link_libraries(${PROJECT_NAME} PUBLIC armadillo ensmallen)
endif() endif()
target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include) target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include)
endif() endif()
target_include_directories(${PROJECT_NAME}
PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph
PUBLIC ${vcglib_devel_SOURCE_DIR}
PUBLIC ${threed-beam-fea_SOURCE_DIR}
# PUBLIC ${MySourcesFiles}
)
if(${MYSOURCES_STATIC_LINK} AND NOT ${USE_POLYSCOPE})
message("Linking statically")
target_link_libraries(${PROJECT_NAME} PRIVATE -static lapack blas gfortran quadmath pthread)
target_link_libraries(${PROJECT_NAME} PUBLIC -static Eigen3::Eigen #[[dlib::dlib]] ThreedBeamFEA #[[${TBB_BINARY_DIR}/libtbb_static.a]])
else()
if(${USE_POLYSCOPE})
target_link_libraries(${PROJECT_NAME} PRIVATE lapack blas gfortran quadmath pthread)
target_link_libraries(${PROJECT_NAME} PUBLIC polyscope xcb Xau)
target_include_directories(${PROJECT_NAME}
PUBLIC ${POLYSCOPE_SOURCE_DIR}
)
endif()
target_link_libraries(${PROJECT_NAME} PRIVATE pthread)
target_link_libraries(${PROJECT_NAME} PUBLIC Eigen3::Eigen #[[dlib::dlib]] ThreedBeamFEA matplot)
endif()
target_link_directories(${PROJECT_NAME} PRIVATE ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
#set(USE_MATPLOT FALSE)
#if(USE_MATPLOT)
##matplot++ lib
download_project(PROJ MATPLOTPLUSPLUS
GIT_REPOSITORY https://github.com/alandefreitas/matplotplusplus
GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${MATPLOTPLUSPLUS_SOURCE_DIR} ${MATPLOTPLUSPLUS_BINARY_DIR})
# add_compile_definitions(MATPLOT_DEFINED)
if(${MYSOURCES_STATIC_LINK})
target_link_libraries(${PROJECT_NAME} PUBLIC "-static" matplot)
else()
target_link_libraries(${PROJECT_NAME} PUBLIC matplot)
endif()
#endif()
download_project(PROJ GSL
GIT_REPOSITORY https://github.com/microsoft/GSL.git
GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${GSL_SOURCE_DIR} ${GSL_BINARY_DIR})
target_link_libraries(${PROJECT_NAME} PUBLIC GSL)
target_include_directories(${PROJECT_NAME} PUBLIC GSL)

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

36
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;
@ -2007,10 +2009,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

12
reducedmodeloptimizer_structs.hpp
View File

@ -222,9 +222,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 +606,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 +615,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 +624,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();

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;

8
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)